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Abstract 

Two  approaches  to  determining  the  ocean  sound  speed  profile  using  measured  acoustic  modal 
eigenvalues  are  examined.  Both  methods  use  measured  eigenvalues  and  mode  dependent 
assumed  values  of  the  WKB  phase  integral  as  input  data  and  use  the  WKB  phase  integral 
as  a  starting  point  for  relating  the  index  of  refraction  to  depth.  Inversion  method  one  is 
restricted  to  monotonic  or  symmetric  sound  speed  profiles  and  requires  a  measurement  of  the 
sound  speed  at  one  depth  to  convert  the  index  of  refraction  profile  to  a  sound  speed  profile. 
Inversion  method  two  assumes  that  the  sound  speed  at  the  surface  and  the  minimum  sound 
speed  in  the  profile  are  known  and  is  applicable  to  monotonic  profiles  and  to  general  single 
duct  sound  speed  profiles.  For  asymmetric  profiles,  inversion  method  two  gives  the  depth 
difference  between  two  points  of  equal  sound  speed  in  the  portion  of  the  profile  having  two 
turning  points,  and  in  the  remainder  of  the  profile  it  gives  sound  speed  versus  depth  directly. 
A  numerical  implementation  of  the  methods  is  demonstrated  using  idealized  ocean  sound 
speed  profiles  numerical  experiments  used  to  test  the  performance  of  the  inversions  using 
noisy  data.  The  two  methods  are  used  to  determine  the  sediment  sound  speed  profiles  in 
two  shallow  water  waveguide  models,  and  inversion  method  one  is  used  to  find  the  sediment 
sound  speed  profile  using  data  from  an  experiment  performed  in  the  Gulf  of  Mexico. 
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Chapter  1 


Introduction 


As  sound  propagates  through  a  medium,  the  sound  field  is  affected  by  the  medium  through 
which  the  sound  is  passing,  and  by  the  nature  and  physical  properties  of  the  medium’s 
boundaries.  This  relationship  between  the  sound  field  and  the  medium  may  be  used  to  predict 
the  properties  of  the  sound  field  as  the  sound  propagates  through  a  medium  with  known 
properties  (the  forward  problem)  or  to  infer  the  properties  of  a  medium  from  measurements 
of  the  sound  field  after  propagation  through  the  medium  (the  inverse  problem).  Using  sound 
we  may  thus  able  to  measure  the  properties  of  media  which  might  otherwise  be  difficult 
to  probe.  In  this  thesis  we  present  a  nonperturbative  inverse  method  based  on  the  WKB 
approximation  to  the  sound  field  which  uses  the  values  of  the  modal  eigenvalues  as  input 
data. 


1.1  Background 

In  ocean  acoustic  tomography,  linear  inverse  techniques  are  applied  to  measurements  of  per¬ 
turbations  in  the  sound  travel  time  which  result  from  sound  speed  variations  along  the  signal 
path.  These  sound  speed  variations  are  then  related  to  the  oceanic  processes  causing  the 
sound  speed  variations.  Acoustic  tomography  techniques  have  been  used  to  study  ocean  fea¬ 
tures  such  as  mesoscale  eddies  (12],  internal  waves  [17],  and  barotropic  motions  [9].  While 
these  processes  may  be  studied  in  other  ways,  tomography  has  the  advantage  of  being  able 
to  provide  frequent  measurements  of  the  ocean  properties  over  large  areas.  Equivalent  mea- 
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surements  made  by  other  means  would  be  prohibitively  expensive.  In  these  tomographic 
studies,  measurements  are  made  using  sound  travelling  over  paths  which  do  not  involve  in¬ 
teraction  with  either  the  ocean  surface  or  bottom,  however;  sound  which  does  interact  with 
the  boundaries  can  be  used  to  study  the  boundaries.  James  Miller  in  his  ScD  thesis  [34], 
for  example,  uses  tomographic  techniques  to  obtain  estimates  of  the  sea  surface  spectra,  and 
acoustic  signals  generated  in  the  water  column  have  been  used  extensively  in  marine  seismol¬ 
ogy  to  study  crustal  structure  and  properties  [16].  For  studying  sound  propagation  in  the 
ocean,  the  properties  of  the  ocean  bottom  are  needed  on  a  much  finer  scale  than  provided 
by  marine  seismological  measurements,  and  the  region  of  interest,  rather  than  extending  to 
several  kilometers  depth,  is  limited  to  the  top  few  hundred  meters  of  the  sediments  at  most 
frequencies  of  interest.  Frisk  et  al.  [19]  use  amplitude  versus  range  data  obtained  using  a 
deep-towed  pulsed  CW  source  and  two  receivers  moored  near  the  bottom  to  infer  geoacoustic 
models  of  the  upper  sediment  layers  in  deep  water.  The  model  is  derived  from  the  measured 
data  using  an  iterative  forward  modelling  technique.  This  method  relies  on  time  gating  to 
separate  the  bottom  reflected  signal  from  the  surface  reflected  signal.  In  shallow  water,  time 
gating  method  is  not  practical  because  the  time  differences  between  signals  reflected  from 
the  surface  and  signals  reflected  from  the  bottom  are  too  small.  Instead,  a  horizontal  array  is 
used  to  measure  the  steady  state  sound  field  which  is  then  decomposed  into  its  eigenfunctions. 
Using  the  eigenvalues  of  the  propagating  sound  field,  Rajan  et  al.  [37]  obtain  the  sediment 
properties  using  a  perturbative  approach. 

Rather  than  using  a  perturbative  technique,  we  use  a  WKB  phase  integral  based  inverse 
method  to  obtain  the  sound  speed  profile  using  as  input  data  the  eigenvalues  of  the  sound  field 
in  the  ocean  acoustic  waveguide.  The  method  may  be  used  to  find  the  sound  speed  profile  in 
the  water  column  and  in  the  upper  layers  of  sediments.  Although  the  sediments  have  enough 
rigidity  to  transmit  shear  waves,  and  we  expect  compressions!  to  shear  wave  conversion  to  take 
place  at  layer  interfaces,  Fryer  [23]  indicates  that  effects  of  compressional  to  shear  conversion 
resulting  from  shear  speed  gradients  within  the  sediment  column  is  small  above  20  Hz.  Since 
we  are  interested  in  the  upper  sediment  layers  where  shear  wave  speed  and  shear  wave  speed 
gradients  are  small,  we  will  treat  the  sediments  as  a  fluid  extension  of  the  water  column  and 
neglect  all  shear  effects.  This  will  enable  us  to  determine  the  compressional  sound  speed 
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profile  in  the  sediments  in  that  portion  of  the  sediments  affecting  sound  propagation  at  the 
frequencies  of  interest  (the  standard  frequency  used  in  the  examples  will  be  220  Hz).  These 
methods  will  not  be  able  to  distinguish  the  low  velocity  zone  that  commonly  occurs  at  the 
water/sediment  interface  for  bottom  materials  such  as  pelagic  clay  [25].  This  inability  to  find 
low  velocity  regions  results  because  the  inversions  are  determining  ray  turning  depths  and 
rely  on  sound  energy  being  turned  by  an  increasing  sound  speed  within  a  region  in  order  to 
find  the  sound  speed/depth  relationship  in  that  region.  The  sound  does  not  turn  within  a 
low  velocity  zone.  Using  an  approach  suggested  by  Vassell  to  study  the  index  of  refraction 
in  optical  waveguides  [43],  we  relate  the  normal  mode  eigenvalues  to  the  sound  speed  profile 
via  the  WKB  phase  integral 

*o  J  [«2(z)  -  stn20„]  /  dz.  (1.1) 

Here  k0  is  the  wavenumber  given  by  u/co  where  u  is  the  frequency  and  Co  is  the  reference 
sound  speed,  6n  is  the  local  propagation  angle  corresponding  to  the  n  th  mode,  and  z\  and 
Z2  are  the  turning  depths  for  the  mode.  The  phase  integral  will  be  treated  as  a  function  of 
the  new  variable  £  =  sin26n.  The  value  of  this  integral  may  be  estimated  based  on  the  mode 
number  and  the  type  of  boundary  interactions  experienced  by  the  mode.  For  example,  the 
phase  integral  for  the  third  mode  will  assume  a  different  value  in  an  environment  giving  two 
turning  points  than  it  will  in  an  environment  where  the  mode  interacts  with  the  surface  and 
has  only  a  lower  turning  point. 

There  are  two  methods  for  obtaining  the  index  of  refraction  as  a  function  of  depth.  The 
first  method  (inversion  method  1)  follows  Vassell’s  treatment  and  can  be  used  for  sound  speed 
profiles  that  increase  monotonically  with  depth  or  are  symmetric  about  a  depth  of  minimum 
sound  speed.  In  this  method,  the  WKB  phase  integral  equation  is  differentiated  to  arrive  at 
am  Abel  integral  equation  which  can  be  used  to  relate  depth  to  index  of  refraction.  From 
knowledge  of  the  sound  speed  at  one  depth,  we  can  determine  the  sound  speed  profile.  If  the 
sound  speed  profile  is  symmetric  rather  than  monotonic,  the  inversion  gives  the  sound  speed 
profile  below  the  depth  of  symmetry,  and  the  symmetry  of  the  profile  is  used  to  generate  the 
upper  portion  of  the  profile. 

The  second  method  (inversion  method  2)  is  based  on  a  method  used  in  quantum  me¬ 
chanical  scattering  for  determining  molecular  potentials  and  is  referred  to  in  the  quantum 
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mechanical  context  as  the  JWKB  or  semi-classical  approximation  [10], [44].  Now  the  WKB 
phase  integral  is  transformed  to  the  form 

itc-jOM*  <‘-2> 

and  /(£)  (referred  to  as  the  inclusion)  is  differentiated  to  produce  the  excursion  which  is  the 
difference  in  depth  between  two  points  of  the  sound  speed  profile  having  the  same  sound 
speed.  In  general  there  is  insufficient  information  to  translate  this  to  a  sound  speed  profile; 
however,  if  the  profile  is  monotonic  or  symmetric,  this  method  gives  the  sound  speed  profile 
directly.  For  an  asymmetric  profile  in  which  some  modes  have  two  turning  points  and  some 
have  one  turning  point,  inversion  method  2  will  give  a  sound  speed  profile  for  the  range  of 
sound  speeds  having  a  single  turning  depth  since  the  surface  is  essentially  the  upper  turning 
depth. 

Munk  and  Wunsch  [35]  derive  an  Abel  transform  based  inversion  method  along  the  lines 
of  our  inversion  method  1;  however,  they  deal  strictly  with  applications  to  the  water  column, 
and  they  approach  the  asymmetric  case  in  a  different  fashion.  Much  of  the  previous  work  in 
inverse  methods  for  finding  the  sound  speed  profile  in  the  ocean  and  the  sediments  has  used 
perturbation  methods  largely  because  the  techniques  of  linear  perturbation  theory  provide  a 
means  to  compute  estimates  of  errors  and  depth  resolution  [33].  With  the  methods  discussed 
here,  we  can  only  give  qualitative  arguments  concerning  resolving  power  and  numerical  ex¬ 
amples  illustrating  accuracy  of  the  methods  and  their  performance  in  a  noisy  environment. 
The  main  advantage  with  our  methods  is  that  they  do  not  require  an  initial  estimate  of  the 
sound  speed  profile. 

1.2  Overview 

Chapter  2  contains  a  review  of  the  essentials  of  normal  mode  theory,  ray  theory,  and  WKB 
theory  in  the  context  of  underwater  acoustics.  Since  the  inversion  methods  presented  here 
are  based  on  the  normal  mode  eigenvalues  as  input  data  and  use  the  WKB  phase  integral 
as  a  starting  point  for  relating  the  index  of  refraction  to  the  eigenvalues,  an  understanding 
of  normal  modes  and  the  WKB  approximation  is  important.  Ray  theory  provides  a  means 
for  building  a  conceptual  picture  of  the  inversion  process  through  the  ray  theory  view  of  the 
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turning  point. 

Chapter  3  contains  the  bulk  of  the  thesis.  We  start  by  grouping  normal  modes  according 
to  the  boundary  interactions  and  the  number  of  turning  points  (we  do  not  allow  for  multiple 
duct  profiles),  and  then  evaluating  the  WKB  phase  integral  for  each  of  the  four  types  of 
mode.  These  phase  integral  values  and  the  associated  eigenvalues  provide  the  input  data  used 
in  generating  a  profile  dependent  functional  relationship  for  the  phase  integral.  The  Abel 
integral  equation  based  inversion  relations  for  the  two  inversion  methods  are  then  derived. 
Following  a  brief  description  of  the  spline  methods  used  in  the  computations,  we  discuss 
the  numerical  implementation  of  the  two  inversion  methods.  Both  inversion  methods  are 
applied  to  the  n2(z)  linear  profile,  the  parabolic  profile,  and  a  bilinear  profile  to  illustrate 
their  capabilities  under  ideal  conditions.  Using  the  n2(z)  profile,  we  then  demonstrate  the 
performance  of  the  inversion  methods  using  data  contaminated  with  random  noise. 

In  chapter  4  we  apply  inversion  method  1  to  the  problem  of  determining  the  sediment 
sound  speed  profile  in  a  shallow  water  waveguide  by  treating  the  sediment  as  a  fluid  extension 
of  the  water  column.  Two  synthetic  models  are  used  to  test  the  inversions:  (1)  a  waveguide 
with  a  terrigenous  bottom,  and  (2)  a  waveguide  with  a  fine  sand  bottom.  We  then  apply  the 
inversion  method  to  data  obtained  in  an  experiment  performed  in  the  Gulf  of  Mexico. 

Chapter  5  summarizes  the  results  of  the  work  and  presents  ideas  for  future  work. 
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Chapter  2 


Normal  Mode  Propagation  and  the 
WKB  Approximation 


2.1  Modal  Propagation 

The  input  data  for  the  inverse  method  considered  in  this  work  are  the  measured  horizontal 
wavenumbers  or  eigenvalues  for  the  modal  description  of  sound  propagation  in  the  medium, 
and  a  brief  review  of  normal  mode  propagation  is  given  here  to  emphasize  the  normal  mode 
ideas  important  to  the  inversion.  Normal  mode  theory  provides  a  Tull  wave’  solution  to 
the  wave  equation  in  the  ocean  waveguide  which  takes  into  account  the  properties  of  the 
medium  in  the  waveguide,  frequency  effects,  and  the  effects  of  the  properties  of  the  waveguide 
boundaries.  The  degree  to  which  the  solution  to  the  wave  equation  provided  by  modal 
theory  realistically  describes  acoustic  propagation  in  the  ocean  depends  on  the  accuracy  to 
which  the  medium  properties  and  boundary  conditions  are  modeled.  We  will  use  models  of 
varying  sophistication  to  illustrate  various  aspects  of  the  normal  mode  description  of  sound 
propagation  pertinent  to  the  inverse  method. 

We  model  the  ocean  as  a  waveguide  consisting  of  a  water  column  bounded  above  by 
the  air-ocean  interface  where  the  acoustic  pressure  is  zero,  and  below  by  a  sediment  bottom 
which  will  be  considered  a  semi-infinite  half-space.  The  sound  speed  is  c(z)  and  the  density 
p  are  taken  to  be  piecewise  constant  throughout  most  of  the  work.  The  following  simplifying 
assumptions  are  used  throughout: 


15 


1.  The  source  is  a  monochromatic  point  source. 

2.  Cylindrical  symmetry  about  the  vertical  axis  through  the  source. 

3.  The  medium  properties  vary  only  in  depth,  i.e.  the  problem  is  range  independent. 
Included  here  is  the  assumption  that  the  water  depth  is  constant. 

4.  No  scattering  from  rough  surfaces  or  inhomogeneities  in  the  medium. 

5.  The  bottom  is  a  fluid  which  may  effectively  be  treated  as  an  extension  of  the  water 
column. 


With  the  assumption  of  cylindrical  symmetry,  the  spatial  part  of  the  acoustic  pressure 
p(r,z,zo),  due  to  a  point  source  located  at  r  =  0  and  z  =  zq  with  harmonic  time  dependence 
exp(-iut),  satisfies  the  inhomogeneous  Helmholtz  equation: 


p(r,z,z0) 


-2 


r  Hr) 


6(z  -  zQ) 


(2.1) 


where  u  is  the  frequency  and  k(z)  =  u>/c(z). 

To  complete  the  problem,  appropriate  boundary  conditions  must  also  be  specified.  At 
any  interface  where  there  is  a  sudden  change  in  the  medium  properties  (the  water-sediment 
interface  for  example)  the  boundary  conditions  of  continuity  of  pressure  and  continuity  of 
normal  particle  velocity  must  be  satisfied.  In  the  range  equation  the  radiation  condition  that 
there  are  no  waves  propagating  inward  from  infinity  will  be  used. 


2.1.1  The  Rigid  Bottom 

We  start  by  considering  the  simplest  model  of  the  ocean  waveguide  consisting  of  a  layer  of 
water  with  sound  speed  c(z),  which  may  be  a  function  of  depth,  and  constant  density,  Po- 
The  layer  is  bounded  above  by  a  pressure  release-surface  (the  air-water  interface)  and  below 
by  a  rigid  bottom.  The  symmetry  and  range  independence  of  the  problem  lead  naturally  to 
use  of  separation  of  variables  in  solving  the  problem.  Letting 

p(r,z,zo)  =  <f>(z)rj>(r),  (2.2) 
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and  substituting  this  into  the  homogeneous  Helmholtz  equation  leads  to  the  separated  equa¬ 
tions 


(2.3) 

and 

(2.4) 

where  -k?  is  a  separation  constant,  and 

= *»(*)  -  e 

(2.5) 

The  total  acoustic  pressure  vanishes  everywhere  on  the  pressure-release  surface  and  a 
sound  wave  incident  on  the  surface  is  completely  reflected  with  a  plane  wave  undergoing  a  ir 
phase  shift.  The  rigid  boundary  condition  is  characterized  by  the  condition  that  the  normal 
particle  velocity  vanishes  on  the  surface  resulting  in  complete  reflection.  For  an  incident  plane 
wave  there  will  be  no  phase  shift  upon  reflection  from  the  rigid  boundary.  In  this  model,  the 
depth  function  4>(z)  satisfies  the  boundary  conditions 


#0)  =  0 


(2.6) 


for  the  pressure-release  surface  and 

d<Kh)  __ 

dz 


(2.7) 


for  the  rigid  bottom. 

Equation  (2.3)  only  has  solutions  for  certain  values  of  the  parameter  7  (the  eigenvalues). 
Because  this  is  a  Sturm- Liouville  problem,  the  set  of  eigenfunctions  <f>i(z),  i  =  1,2,3,..., 
forms  a  complete  orthonormal  set  of  square  integrable  functions  on  the  interval  0  <  z  <  h 
implying  that  [45] 

fk  MZ)MZ) 

Jo  Po 

where  the  l//>o  weighting  is  included  in  analogy  to  the  results  in  problems  where  the  density 
is  a  function  of  depth. 

Up  to  this  point  we  have  allowed  the  sound  speed  to  be  any  function  of  the  depth  z.  For 
definiteness  in  illustrating  the  solution  of  the  problem,  we  now  specify  c(z)  =  constant.  The 


(2.8) 
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boundary  conditions  lead  to  the  equation 


7  n  ~ 


(2n  -  1)t 
2  h 


for  the  vertical  wavenumber  7„.  Using  a  density  po  =  1  and  normalizing  gives  for  the 
eigenfunctions 

<M*)  =  y/fsin  [(*2  -  *r2J1/2*]  (2.9) 

where  the  horizontal  wavenumbers  kTn  are  given  by 

r(2n-  1)*T2 


I k2  =fc2- 


2h 


(2.10) 


Expanding  the  solution  of  equation  (  2.1)  in  terms  of  the  eigenfunctions 

OO 

p(r,z,z0)  =  £  </>n(z)i/>n(r), 


(2.11) 


n=l 


substituting  expansion  into  equation  ( 2.1),  and  using  the  orthonormality  of  the  depth  func¬ 
tions  gives 

/  r)  1  //«/i  o 

(2.12) 


^tMr)  l<fV>n(r)  L  u2 
dr2  r  dr 


+  — +  *?„tfw(r)  =  -h(r)Mzo) 


which  is  a  Bessel  equation  of  order  zero  for  which  the  solution  is  a  linear  combination  of  Hankel 
functions.  We  use  the  radiation  condition  that  there  are  no  waves  propagating  inward  from 
infinity  to  eliminate  the  term  involving  Hq2\  and  the  solution  to  the  radial  equation  is 


xl>n(r)  =  iff<M*o)#o1)(fcr„T'). 


(2.13) 


The  solution  to  equation  ( 2.1)  is  then  expressed  as  a  sum  of  the  normal  modes 

f)  *  OO 

p(r,z,zo)  -  -P^sin(7nz)sin(7r,zo)^o1)(fcr„r). 


(2.14) 


At  distances  of  more  than  a  few  wavelengths  from  the  source,  we  can  use  the  asymptotic 
approximation  to  the  Hankel  function  [2] 


~  \l exP[i(kTr  -  */4)] 


(2.15) 


The  normal  mode  solution  of  equation  (  2.1)  given  by  equation  (  2.14)  describes  the 
sound  field  in  terms  of  a  sum  of  travelling  cylindrical  waves  propagating  outward  from  the 
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origin.  Each  of  these  modes  is  a  standing  wave  in  the  vertical  direction  with  a  distinct  depth 
distribution  of  pressure  described  by  the  eigenfunction  4>{z)  and  the  vertical  wavenumber  yn. 
Although  the  summation  of  equation  ( 2.14)  is  infinite,  for  a  sufficiently  large  mode  number  N , 
the  horizontal  wavenumber  will  become  imaginary,  and  the  solution  will  rapidly  decay  with 
range.  For  a  receiver  located  sufficiently  far  from  the  source,  the  sound  field  will  be  adequately 
described  by  a  finite  sum  of  propagating  modes.  It  is  these  horizontal  wavenumbers  from 
this  description  of  the  sound  field  that  will  provide  the  input  information  for  the  inversion 
method  described  in  the  next  chapter. 


2.1.2  Normal  modes  in  a  waveguide  with  a  penetrable  bottom 


If  the  bottom  of  the  waveguide  is  not  a  rigid  reflector,  some  of  the  sound  energy  will  penetrate 
into  the  bottom,  and  the  solution  to  the  Helmholtz  equation  must  be  found  both  in  the  water 
and  in  the  bottom.  The  approach  taken  here  is  different  from  the  previous  section  in  that 
the  Hankel  transform  is  used  to  obtain  the  solution;  however,  the  essential  ideas  related  to 
the  propagating  modes  remain  the  same. 

Starting  once  again  with  the  inhomogeneous  Helmholtz  equation  (equation  (  2.1)),  we 
introduce  the  zero  order  Hankel  transform  pair 


p(r,  z,zo)  = 

roo 

/  p{kr,z,zo)Mkrr)krdr 

Jo 

(2.16) 

p(kr,Z,Zo) 

[°° 

=  /  p(r,z,ZQ)J0{krT)rdr 

Jo 

(2.17) 

where  Jo(krr)  is  the  Bessel  function  of  order  zero. 

Multiplying  equation  (2.1)  by  Jo(kTr)rdr  and  integrating  from  0  to  oo  (see  [3]  for  details 
of  the  integration),  we  find  that  the  transform  variable  (or  depth  dependent  Green’s  function) 
must  satisfy  [20] 

+  fc2(z)  -  kij  p(kr,z,zo)  =  -2 6(z  -  *o).  (2.18) 

Let  pi  and  fo  be  linearly  independent  solutions  of  equation  (2.18)  chosen  such  that  pi  satisfies 
the  surface  boundary  condition,  and  p?  satisfies  the  bottom  boundary  condition.  Then  the 
transform  function  p{kr,z)  is  given  by  [20] 


p(kr,z,zo)  = 


-2pi{kr,z)p2{kT,z0) 

W(zo) 


0  <  z  <  Zo 


(2.19) 
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Zo  <  z 


(2.20) 


where  the  prime  indicates  differentiation  with  respect  to  z. 

Next  the  Bessel  function  is  expressed  in  terms  of  Hankel  functions  as 

J0(krr)  =  i  [H^(krr)  +  42)(krr)]  (2.22) 

and  use  the  relation 

H^ie-^krr)  =  -H^](krr)  (2.23) 

to  allow  the  range  of  integration  for  the  transform  integrals  to  be  extended  over  the  whole 
real  kr  axis  from  -oo  to  oo  when  evaluating  the  transform  integral  for  p(r,z).  The  solution  is 
then  found  by  contour  integration  methods  [3],  and  may  be  expressed  as  the  sum  of  a  discrete 
and  a  continuous  portion  given  by  [20] 

P(r,z )  =  i*YlMzc)<t>n(*)41\krnr)  +  /(r)  (2.24) 

where  the  eigenfunctions  and  eigenvalues  satisfy 

+  *!«  -  *?.)  *.(*>  = 0 

along  with  the  surface  and  bottom  boundary  conditions. 

The  part  of  the  solution  representing  the  continuum  (/(r)  or  the  branchline  contribution) 
is  rapidly  attenuated  with  range  and  will  be  neglected,  since  we  are  only  concerned  with 
measurements  made  more  than  a  few  wavelengths  from  the  source.  We  note  that  the  dis¬ 
crete  portion  of  the  solution  is  of  the  same  form  as  for  the  rigid  bottom  case,  and  that  the 
eigenvalues  have  the  same  interpretation,  although  they  will  take  on  different  values  since 
the  sound  is  propagating  in  both  media,  and  the  characteristic  equation  will  be  different  .In 
general  the  lowest  medium  will  be  taken  as  a  half  space  extending  to  infinity  and  a  radiation 
condition  will  be  imposed  resulting  in  a  solution  that  decays  with  depth  in  the  halfspace. 
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2.1.3  Analytically  solvable  profiles 

We  next  consider  two  specific  sound  speed  profiles  in  order  to  illustrate  the  concept  of  turning 
depth.  A  turning  depth  is  a  depth  at  which  the  vertical  wavenumber  vanishes,  leading  to 
a  fundamental  change  in  the  behaviour  of  the  solution  to  the  Helmholtz  equation.  In  the 
region  where  the  vertical  wavenumber  is  greater  than  zero,  the  solution  is  oscillatory,  while 
in  the  region  for  which  the  vertical  wavenumber  is  less  than  zero,  the  solution  will  decay 
with  depth.  Depending  on  the  sound  speed  structure,  these  turning  points  may  be  located 
in  either  the  water  column  or  in  the  sediment.  The  sound  speed  profiles  we  will  consider  are: 
(1)  the  square  of  the  index  of  refraction  a  linear  function  of  depth,  and  (2)  the  square  of  the 
index  of  refraction  parabolic.  In  both  these  cases  solution  of  the  depth  equation  gives  closed 
form  expressions  for  the  eigenvalues.  Additionally,  we  will  later  find  that  these  profiles  may 
be  inverted  analytically  by  the  methods  presented  here,  and  consequently  they  will  be  the 
standard  examples  for  illustrating  results  and  verifying  the  numerical  methods. 

The  n2(z)  linear  waveguide 

For  simplicity  we  consider  a  semi-infinite  ocean  with  a  pressure-release  surface  and  an  index 
of  refraction  which  satisfies  the  relation 

n2(z)  =  1  -  az.  (2.26) 

In  lieu  of  a  bottom  boundary  condition,  we  apply  a  radiation  in  depth  requiring  that  the 
solution  remain  finite  at  infinity. 

Proceeding  as  for  the  penetrable  bottom  case,  the  homogeneous  part  of  equation  ( 2.18) 
becomes 

^  +  [*o(l  ~  az)  -  p(k r,  z,  zo)  =  0  (2.27) 

Making  the  change  of  variables  [8] 

H  =  (a*2)-1'3;  to  =  H*(k?  -  *2);  t  =  t0  +  z/H ,  (2.28) 

equation  ( 2.27)  is  transformed  into  the  Airy  equation 

=  tp{t)  (2.29) 
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which  has  as  solutions  the  Airy  functions  Ai(t)  and  Bi(t)  [2].  (Here  functional  dependence 
other  than  dependence  on  t  has  been  surpressed.)  We  eliminate  Bi(t)  in  order  to  satisfy  the 
radiation  condition,  and  the  eigenvalues  are  found  using  the  pressure  release  condition  which 
requires  At(to)  =  0  from  which  we  see  that  to  =  -yn  where  the  yn  are  the  zeros  of  the  Airy 
function.  With  tn  =  z/ H  -  yn  and  ton  =  zo/H  -  yn  the  sound  field  is 

iff  ^  Ai(tn)Ai(ton)Ho\krr)  ,n 

'Kr'2)=«? — — '  ,2'!0) 

The  eigenvalues  are  given  by 

k2n  =  k2-yJH2.  (2.31) 


The  parabolic  profile 

In  considering  the  parabolic  profile,  we  take  the  ocean  as  infinite  in  the  vertical  direction 
so  that  there  are  no  boundaries,  and  the  radiation  condition  is  used  to  give  a  solution  that 
remains  finite  at  ±oo.  The  index  of  refraction  satisfies 

n2(z)  =  1  —  a2z2,  (2.32) 

and  the  homogeneous  part  of  equation  ( 2.18)  is 

+  \ko(l  -  °2*2)  -  kr]  Pikr,z,zo)  =  0  (2.33) 

The  solution  to  this  are  given  in  terms  of  Hermite  polynomials  as  [42] 

P(kr ,  z ,  zq  )  =  (2  nn\)-'l2{akol*)xlAe-^Hn(y/^>z)  (2.34) 

where  the  Hn{-)  are  the  Hermite  polynomials.  The  eigenvalue  condition  gives 

fcn  =  Ml-<»(2n  +  l)/*o)1/2;  n  ss  0, 1,2 .  (2.35) 


2.2  Rays 

Ray  theory  is  a  geometric  solution  of  the  wave  equation  which  is  correct  for  high  frequencies 
and  provides  a  simple  physical  description  of  the  solution  in  terms  of  the  paths  along  which 
the  acoustic  energy  is  refracted.  The  ray  approximation  is  widely  used  for  studying  sound 
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propagation  in  underwater  acoustics,  and  provides  a  convenient  method  for  visualizing  the 
propagation  paths  of  the  sound  energy  and  the  depth  of  penetration  of  the  sound  energy  into 
the  bottom.  We  will  briefly  illustrate  the  connection  between  the  normal  mode  view  and  the 
ray  picture  to  the  extent  that  it  helps  in  understanding  the  inversion  technique.  A  number 
of  references  are  available  which  provide  detailed  treatments  of  ray  theory  in  underwater 
acoustics  [7], [6]. 

2.2.1  The  Equations  of  Ray  Theory 

The  equations  for  ray  acoustics  [7]  are  obtained  from  the  Helmholtz  equation  in  vector  form 
using  pressure  as  the  variable 

V2p  +  k\x)p  =  0,  (2.36) 

where  x  is  the  position  vector  and  k(x )  =  u/c(x),  by  tahing  the  solution  to  be 

p{x)  =  A{x)exp{ikoS{x)).  (2.37) 

Here  k0  =  u>/co  where  cq  is  a  reference  sound  speed  such  as  the  minimum  sound  speed  in  the 
waveguide  (or  the  sound  speed  at  the  source),  A(x)  is  the  amplitude  of  the  sound  waves,  and 
k0S(x)  is  the  phase,  with  the  function  S( x)  referred  to  as  the  eikonal. 

Substituting  the  assumed  solution  (  equation  (2.37))  into  the  Helmholtz  equation  gives 

V2A  +  »fc0(2VA  •  VS  +  AV2S)  +  klA  [n2(x)  -  (VS)2]  =  0.  (2.38) 

The  equations  of  ray  theory  are  obtained  from  equation  ( 2.38)  in  the  limit  as  ko  — ♦  oo  by 
equating  the  real  and  imaginary  terms  to  0,  and  neglecting  the  V2A  compared  to  the  real 
term  containing  k%A  (that  is  we  assume  V2A/A  <  k%).  The  resulting  equations  are  the 
eikonal  equation 

(VS)2  =  n2(x),  (2.39) 

and  the  transport  equation 

2 VA  •  VS  +  AVS  =  0.  (2.40) 

The  eikonal  equation  defines  the  geometry  of  the  rays  which  are  lines  perpendicular  to  the 
surfaces  of  constant  phase,  that  is  S(x)  =  constant  .  These  ray  trajectories  can  be  computed 
to  determine  the  path  of  the  sound  energy. 
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2.2.2  The  Connection  Between  Rays  and  Modes 

C.T.  Tindle  and  K.M.  Guthrie  [24]  [41]  have  illustrated  the  connection  between  normal  mode 
theory  and  the  ray  theory  approximation.  Considering  a  deep  water  sound  speed  profile  with 
a  minimum  sound  speed  at  some  depth  below  the  surface,  then  for  a  given  mode  number  n, 
a  ray  with  the  same  turning  point  depths  as  mode  n  can  be  defined  by 


Cq  _  u>_ 
sin  0O  kTn 


(2-41) 


where  the  angle  $o  is  the  angle  of  incidence  for  the  ray  as  it  crosses  the  channel  axis  (i.e.  the 
sound  speed  minimum),  and  kTn  is  the  horizontal  wavenumber  for  mode  n.  The  ray  effects 
are  duplicated  by  summing  over  a  group  of  modes,  and  each  group  will  manifest  itself  as  an 
energy  pulse  travelling  along  the  ray  path  which  corresponds  to  the  ray  defined  by  the  angle 
Oq  for  the  central  mode  in  the  group  of  modes.  Thus  we  can  construct  a  ray  corresponding 
to  the  group  of  modes;  however,  this  ray  will  not  necessarily  be  the  same  as  the  physical 
ray  path  along  which  the  sound  energy  actually  travels  between  the  source  and  receiver.  By 
using  this  ’dew  of  rays  as  the  result  of  modal  interference,  we  can  construct  a  picture  of  the 
path  followed  by  the  sound,  and  its  turning  depths  within  the  framework  of  normal  modes. 


2.3  The  WKB  Approximation 

WKB  theory,  like  ray  theory,  provides  an  approximate  solution  to  the  wave  equation  in 
the  high  frequency  regime;  however,  the  WKB  method  takes  frequency  information  into 
account  in  the  amplitude  function  as  well  as  in  the  phase,  and  it  accounts  for  frequency 
dependent  effects  such  as  diffraction  which  are  ignored  in  standard  ray  theory.  In  general,  the 
WKB  method  is  a  means  for  finding  an  approximate  solution  to  a  linear  ordinary  differentia] 
equation  in  which  the  highest  order  derivative  is  multiplied  by  a  small  parameter  [5]  and  its 
use  is  not  limited  to  underwater  acoustics.  In  the  problem  at  hand,  the  differential  equation 
of  interest  is  the  equation  for  the  depth  dependent  Green’s  function  (equation  ( 2.18)) 

~~7T~ Z°^  +  *oh2(*)  -  k?/kl]p(kr,z,zo)  =  0 
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where  ko  =  w/co,  Co  is  the  sound  speed  minimum,  and  n(z)  is  the  index  of  refraction.  For 
later  convenience  we  make  the  following  change  in  variables 

k(z)  =  n2(z),  ( =  k2Jkl  (2.42) 


and 

=  «(*)  -  f. 


(2.43) 


Equation  (  2.3)  then  becomes 

d?p(kT,z,zo) 

d P 


+  koQ2{z)P(kT,z,zo)  =  0. 


This  equation  is  solved  by  expanding  p  as  [22] 


(2.44) 


p(kr,z,zo)  =  exp 


(2.45) 


where  the  lower  limit  zq  is  some  conveniently  chosen  constant  depth.  Substituting  this  ex¬ 
pansion  into  equation  ( 2.44  )  and  setting  the  coefficients  of  successive  powers  of  k0  equal  to 
zero,  we  obtain 


yo(z)  =  ±iQ(z)  (2.46) 

=  -£yMy„_M;  v~  1|2»3,...  (2.47) 

<i=0 

from  which  the  series  in  the  exponent  of  equation  ( 2.45)  may  be  determined.  This  series  will, 
in  general,  be  asymptotic  and  not  convergent.  Retaining  the  first  three  terms  and  using  a 
convenient  normalization  gives  the  second  order  WKB  approximation 


p(fcr, z, zo)  =  (fc0Q(z))-1/2  exp  |±  J  ^1  +  k0Q(z)dz 

where 

«  =  «j«*or3/!  gp  (*o<j(z)r,/! . 

The  steps  leading  to  equation  ( 2.48)  are  justified  if 


(2.48) 


I  1/2*0  2  I  <  I  yi*0  1  t  <  I  i/0  I 


(2.49) 


or 


<o  ]  <1  jz  (koQ(z)rl 


<2. 


(2.50) 
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If 


I  [*  koe0Q(z)dz  |  <  2  (2.51) 

Jzo 

then  we  can  make  the  physical  optics  approximation  [5]  of  retaining  only  the  first  two  terms, 
and  the  two  linearly  independent  solutions  of  equation  ( 2.3)  are 

P(kT,z,z0)  =  kol/2(n(z)  -  t)~1,Aezp  ±  ik0  f  (k(z)  -  £)l/2  dz  .  (2.52) 

Jzo 


The  WKB  method  has  provided  an  approximate  solution  to  equation  (2.18)  which  takes 
into  account  sound  speed  profile  and  frequency  information,  and  which  allows  for  frequency 
dependent  effects  such  as  diffraction.  These  WKB  modal  solutions  can  then  be  used  in  the 
modal  sum  (equation  (2.24))  to  approximate  the  full  wave  solution.  It  is  well  known  that  the 
WKB  solution  fails  at  the  turning  points  i.e.  the  points  at  which  Q(z )  =  («(*)  -  £)  =  0.  As 
noted  earlier,  the  behaviour  of  the  solution  differs  on  the  two  sides  of  the  turning  point  with 
oscillatory  behaviour  above  the  turning  point  and  an  exponentially  decaying  behaviour  below 
the  turning  point.  There  is  an  extensive  literature  concerning  the  problem  of  connecting  the 
solutions  through  the  turning  point  (see  [5]);  however,  we  are  only  interested  in  the  phase 
memory  of  the  solution  i.e.  the  integral  in  the  exponential  for  the  region  above  the  turning 
point,  and  not  in  generating  approximations  to  the  sound  field. 
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Chapter  3 


The  Inversion 


We  approach  the  inversion  by  first  grouping  modes  according  to  the  boundary  interactions 
they  undergo  and  the  number  of  turning  points.  Based  on  these  factors  the  WKB  phase  inte¬ 
gral  is  evaluated  for  each  type  of  mode  resulting  in  an  expression  which  depends  only  on  the 
mode  number  in  the  case  of  modes  which  do  not  interact  with  the  bottom,  and  on  mode  num¬ 
ber  and  the  phase  of  the  bottom  reflection  coefficient  for  bottom  interacting  modes.  Treating 
these  expressions  as  functions  of  the  variable  £,  we  proceed,  for  each  inversion  method,  to 
obtain  an  Abel  integral  equation  which  enables  us  to  extract  a  relationship  between  the  index 
of  refraction  and  either  the  depth  (inversion  method  1)  or  excursion  (inversion  method  2). 
Both  inversion  methods  are  applied  to  three  prototypal  sound  speeed  profiles  to  demonstrate 
their  performance  under  ideal  conditions.  The  performance  of  the  inversions  is  then  tested 
using  one  profile  (the  n2(2)  profile  )  with  random  errors  added  to  the  eigenvalues. 

3.1  Evaluating  the  Phase  Integral 

There  are  four  basic  types  of  normal  modes  depending  on  the  region  where  the  normal  modes 
are  concentrated  i.e.  the  region  where  y/n  -  £  is  real  [7]  (sound  speed  profiles  with  multiple 
sound  ducts  are  excluded  from  this  discussion): 

1.  0  <  2  <  zilavtT.  The  region  is  bounded  above  by  the  pressure  release  surface  and  below 
by  a  turning  point. 
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2.  ztuPPtr  <  2  <  zilom„.  The  region  is  bounded  above  by  a  turning  point  at  Ztupper  and 
below  by  a  turning  point  at  ztloulcr.  On  the  ocean  surface  and  bottom,  the  sound  field 
is  exponentially  small.  This  case  would  include,  for  example,  SOFAR  propagation. 

3.  0  <  z  <  H.  The  region  is  bounded  above  by  the  ocean  surface  and  below  by  the  ocean 
bottom. 

4.  ztupptr  <  z  <  H.  The  region  is  bounded  above  by  a  turning  point  and  below  by  the 
ocean  bottom. 


Because  we  are  interested  in  the  properties  of  the  sediments  as  well  as  the  water  column,  we 
include  in  the  categories  having  lower  turning  points,  those  instances  in  which  the  sound  is 
not  totally  reflected  at  the  interface,  but  turns  within  the  sediment  layers. 

For  convenience,  let  zt,  =  z«louiir.  zti  ~  ztupptT,  and  7  =  (k  -  O1^2-  The  WKB  solution  to 
the  depth  equation  for  the  transform  variable  p(z)  is 


p(z)  =  —  Ci exp  (ikoj*  S^)+  ^exp  (~ik0  jf  1  7^)]  ; 


zta  <  z<  ztl  (3.1) 

p(z)  =  -^L^Cacxp  ^-*0  jT  I  7  I  ;  Ztl  <  z.  (3.2) 

The  two  terms  in  equation  (3.1)  represent  upgoing  and  downgoing  waves  in  the  sound  channel 
while  equation  ( 3.2)  represents  an  exponentially  attenuated  wave  below  the  turning  depth. 
The  constants  are  related  by  [8] 


Ci  =  Caexp(*jr/4);  C2  =  Czexp(-ir  j  A) 


with  the  tt/4  factor  accounting  for  the  phase  associated  with  turning  at  the  turning  depth. 
Our  interest  is  above  the  turning  depth  where  the  solution  is  oscillatory,  and  that  solution 
can  be  written 

p(z)  =  ~  cos  (k0J'  1  7  dz-  *74^  (3.3) 

from  which  the  phase  integral  will  be  evaluated  using  the  boundary  conditions.  We  are  not 
interested  in  the  value  of  the  constants  C3,  and  will  deal  only  with  the  phase  integral. 
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For  the  first  mode  type  (case  1),  the  pressure  release  condition  at  the  ocean  surface 
requires  that  p(0)  =  0.  Setting  the  cosine  term  equal  to  zero  requires 


*o  /  '  Indz  =  (n  -  1/4)jt;  n  =  l,2,3,... 
Jo 


(3.4) 


For  the  modes  with  two  turning  points  (case  2),  we  assume  that  the  turning  points  are  far 
enough  away  from  the  boundaries  that  interactions  with  the  boundaries  are  not  of  concern. 
Expanding  the  cosine  into  the  sum  of  exponentials 


cos 


(ko  J  1  7n  dz-  =  ^exp  (ik0  J  '  yndz-  iit/Aj 

+  ^exp  ( ~ik°Jz  *  7n  dz  +  tx/4^  . 


(3.5) 


The  first  term  is  a  wave  propagating  in  the  positive  z  direction  (down),  while  the  second  term 
is  a  wave  propagating  in  the  negative  z  direction  (up)  i.e.  a  wave  which  was  reflected  at  the 
upper  turning  point.  At  the  turning  depth,  the  wave  is  reflected  with  a  x/2  phase  change. 
Taking  the  ratio  of  the  two  terms  and  setting  the  ratio  equal  to  exp[i7r(2m  +1)]  gives 


*0  /  1  In  dz  =  (n  -  1/2)  r;  n  =  l,2,3 .  (3.6) 


with  m  =  (n  -  1)  used  to  provide  indices  starting  at  1. 

For  the  case  3  we  consider  a  mode  interacting  with  the  pressure-  release  surface  above 
and  totally  reflected  at  the  bottom.  Expanding  the  field  into  up  and  downgoing  waves  again 
gives 


p(z)  =  A  exp  (i*o  Jq  7n  +  B  exp  (- ik0  J  7„  dz'j 


(3.7) 


Evaluating  the  terms  at  the  two  boundaries  in  terms  of  the  surface  and  bottom  reflection 
coefficients 


p  _  m  i  _  a 

1  mi  b 

(3.8) 

and 

=  21 /„  T"*)' 

(3.9) 

Eliminating  the  ratio  A/ B  leaves 

/  rH  \ 

Ri,Rtexp\2i  J  "fndzj  =  1. 

(3.10) 
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Letting  R,  =  exp(i2$a)  and  Rb  =  e:cp(  t2$f>),  gives 

[H 

+  $b  +  /  Indz  =  (n  -  1) w;  n  =  1,2,3, ... . 
Jo 

For  a  pressure  release  surface  =  -ir/2  and 


(3.11) 


f 


1 Indz  =  (n-  1/2) r  -  $b 


(3.12) 


where  is  the  phase  of  the  reflection  coefficient  for  total  internal  reflection  at  the  angle  of 
incidence  for  the  given  mode.  In  the  case  of  the  Pekeris  waveguide  (  an  isovelocity  water 
column  over  an  isovelocity  bottom  half-space)  this  is  [11] 


$6 

bi 


bipoco 

arctan - — 

Pl  Cl  COS  UQ 


(3.13) 

(3.14) 


where  the  p  are  the  densities,  c  are  the  sound  speeds,  and  Oq  is  the  angle  of  incidence  in  the 
upper  medium.  The  subscript  o  refers  to  the  upper  medium,  and  the  subscript  l  refers  to  the 
lower  medium. 

In  case  4  the  wave  reflected  at  the  upper  turning  depth  lags  in  phase  by  7t/2  behind  the 
incident  wave  [7],  and  =  -ir/4  giving  for  the  phase  integral 


rH 

ko  Indz  =  (n  -  3/4)t  -  n=  1,2,3,.. 

Jzt 


(3.15) 


3.2  Derivation  of  the  Inversion  Relations 

We  designate  the  WKB  phase  integral  divided  by  the  wavenumber  ko  as  F(£).  The  variable 
f  is  physically  equal  to  sin2  6  where  6  is  the  angle  of  incidence  of  the  ray  corresponding  to 
the  mode  having  the  eigenvalue  koy/%.  Viewing  the  modes  as  a  consequence  of  constructive 
interference  of  rays,  we  treat  £  as  a  continuous  variable  with  the  measured  values  of  the 
horizontal  wavenumber  and  the  WKB  phase  integral  defining  values  of  the  function  F(£)  at 
selected  points.  An  Abel  integral  equation  which  can  be  used  to  relate  depth  to  the  index 
of  refraction  (or  sound  speed  depending  on  the  normalization  used)  may  be  obtained  by  two 
methods.  The  first  approach  is  applicable  to  sound  speed  profiles  which  are  monotonically 
increasing  with  depth  or  symmetric  with  respect  to  the  sound  speed  minimum.  The  second 
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approach  allows  determination  only  of  the  distance  between  points  on  the  two  branches  of 
the  sound  speed  profile  having  the  same  sound  speed  when  the  profile  is  asymmetric.  If  some 
of  the  modes  for  the  asymmetric  profile  have  two  turning  points  and  some  have  one  turning 
point,  then  for  sound  speed  range  for  which  there  are  two  turning  points  the  result  will  be  as 
described,  but  for  the  sound  speed  range  with  one  turning  point,  the  result  will  be  the  sound 
speed  profile.  For  a  sound  speed  profile  which  is  monotonically  increasing  with  depth  or  is 
symmetric  about  the  sound  speed  minimum,  the  second  method,  like  the  first,  provides  depth 
as  a  function  of  index  of  refraction  .  We  will  refer  to  the  first  method  which  is  applicable 
to  monotonic  and  symmetric  profiles  as  inversion  method  1,  and  we  will  refer  to  the  second, 
more  general,  method  as  inversion  method  2. 

3.2.1  Derivation  of  inversion  method  1 

Considering  the  case  of  a  profile  which  is  monotonically  increasing  with  depth,  the  phase 
integral  of  the  WrKB  solution  to  the  wave  equation  (assuming  constant  density)  is  written  as 

(3.16) 

Here  k0  is  ui/cq,  c0  is  the  value  of  the  sound  speed  at  the  surface,  and  fcT„  is  the  horizontal 
wavenumber  for  the  nth  mode.  The  equation  has  been  normalized  by  dividing  through  by  fc0, 
and  the  quantity  which  will  ultimately  result  from  the  inversion  is  the  index  of  refraction  n(z). 
Here  zt  is  used  to  denote  the  single  turning  depth.  If  the  sound  speed  profile  is  symmetric, 
the  phase  integral  is  taken  between  upper  and  lower  turning  points,  and  the  symmetry  of  the 
problem  allows  the  integral  to  be  split  into  two  equal  pieces  with  the  resulting  factor  of  two 
taken  to  the  right-hand  side  of  the  equation. 

With  the  change  of  variables  from  the  previous  chapter 

K(z)  =  n2(z);  (  =  (kr„/k 0)2, 

equation  (3.16)  is  written  as 

(3.17) 

Jo 

Following  the  treatment  of  Vassell  [43],  this  can  be  transformed  into  an  Abel  integral  equation 
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by  differentiating  with  respect  to  £  giving 


dFtf) 


1  /** 

=  ^  =  -2/0 


dz 


n  2  h  w2)-{],/:' 

The  independent  variable  z  is  eliminated  in  favor  of  the  dependent  variable  k  by  letting 


(3.18) 


=  -  £  G(K)dt 


(3.19) 


where  G(k)  is  an  unknown  function,  and  in  fact  we  will  not  have  any  reason  to  determine 
<j(k)  explicitly. 

Since  k(z )  is  the  index  of  refraction  squared,  at  the  reference  depth  k(0)  =  1,  and  at  the 
turning  depth  n(zt )  =  Equation  (3.18)  becomes 

1  f 1  G(K)dK 


jj  ( r\  1  /  G(K>dK 

which  is  an  Abel  integral  equation  and  can  be  inverted  [28]  to  give 

ru ,  2  d  H(m 

Using  equation  (  3.19)  and  integrating  gives  the  relation 

H(Qdt 


(3.20) 


(3.21) 


=  -2/ 

jt  Jk 


(3.22) 


K  -  «),/J 

which  can  be  used  to  calculate  k(z)  from  the  H(£)  obtained  by  differentiating  the  function 
F(t)  defined  by  the  measured  horizontal  wavenumbers  and  associated  values  of  the  WKB 
phase  integral.  The  value  of  the  sound  speed  at  one  depth  is  now  used  to  convert  the  index 
of  refraction  to  a  sound  speed  profile. 


3.2.2  Derivation  of  inversion  method  2 

In  the  case  of  an  asymmetric  sound  channel  with  two  turning  points,  we  are  able  to  obtain  a 
relation  for  the  distance  between  the  points  on  the  two  branches  of  the  sound  speed  profile 
having  the  same  sound  speed.  The  information  contained  in  the  phase  integral  is  insufficient 
to  obtain  the  depths  separately.  An  analagous  problem  in  quantum  mechanics  is  the  deter¬ 
mination  of  a  potential  from  bound  state  information  and  the  technique  we  use  is  referred 


32 


to  as  the  JWKB  semi-  classical  approximation  in  the  quantum  mechanical  context  [10], [44]. 
Define  the  inclusion  as 

/(O  =  f"  (k- ()  dz  (3.23) 

Jn 

and  the  excursion  as 

X(0  =  r  dz~zx-  z2.  (3.24) 

J*  t 

It  is  the  excursion  which  we  will  be  able  to  obtain  from  measurements  of  the  horizontal 
wavenumbers  and  associated  values  of  the  phase  integral.  For  a  monotonic  profile,  the  second 
turning  point  is  at  the  surface  ( z  =  0),  and  the  excursion  will  give  the  index  of  refraction 
versus  depth  directly,  and  in  asymmetric  channels  the  surface  is  the  upper  turning  point  for 
modes  interacting  with  the  surface  so  that  the  excursion  will  be  the  sound  velocity  profile  for 
a  portion  of  the  profile. 

Define  the  adjoint  fractional  integral  as  [39] 

*Wb-' 7(x)  =  jf  (C  -  *r 1  /(C¥C;  Re  v>0.  (3.25) 

Applying  this  to  both  sides  of  the  WKB  phase  integral  equation  with  v  =  -1/2  gives 

^m“(r-o,/2[£2(«-n,/2^]  (3.26) 

Figure  (3-1)  shows  the  region  over  which  the  square  of  the  index  of  refraction  is  being 
integrated.  Interchanging  the  order  of  integration,  it  is  apparent  that  t.he  integral  over  £'  will 
be  zero  for  £  greater  than  k.  This  leaves 


(3.27) 


which  reduces  to 


'«>=/>-  (328) 

Taking  the  derivative  of  both  sides  of  this  equation  with  respect  to  £  results  in  an  expression 


for  the  excursion 


(3.29) 


Thus  from  the  values  of  F(£),  we  determine  the  inclusion  for  the  sound  speed  profile,  and 
use  the  inclusion  to  find  the  excursion. 
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Figure  3-1:  The  integration  over  is  carried  out  over  the  cross-hatched  region  extending 
from  £  to  the  curve  k(z). 

3.2.3  Examples  of  analytically  invertible  profiles 

The  n2(z)  linear  profile  and  the  parabolic  profile  are  two  examples  of  sound  speed  profiles 
which  are  analytically  invertible  with  these  techniques,  and  which  have  well  known  expressions 
for  the  eigenvalues.  These  profiles  will  serve  as  test  cases  for  checking  the  accuracy  of  the 
intermediate  steps  in  the  numerical  inversions. 

The  n2(z)  linear  profile 

With  the  square  of  the  index  of  refraction  given  by 

k(z)  =  n2(z)  =  1  -  oz, 

the  WKB  phase  integral  can  be  evaluated  directly 

m  =  l"  ((1  -  «)  -  ()'/2dz  =  £(1  -  (3.30) 
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Using  the  first  inversion  method  we  take  the  derivative  with  respect  to  £  giving 


-1 


(3.31) 


The  integral  expression  for  depth  in  terms  of  the  index  of  refraction  is 


(3.32) 


which  upon  integrating  and  solving  for  k  returns  the  original  expression  for  the  square  of  the 
index  of  refraction. 


The  parabolic  profile 

For  the  parabolic  profile  the  square  of  the  index  of  refraction  is  given  by 

k(z)  =  n2(z)  =  1  -  a2z2,  (3.33) 

and  the  inclusion  is 

HO  =  I"'  ((1  -  «V)  -  i]iz  =  ±  (1  -  if2 .  (3.34) 

Taking  the  derivative  with  respect  to  £  gives  the  excursion 

“  =  x(0  —  ~(1  -  01/2*  (3-35) 

The  depth  of  the  turning  points  is  given  by  the  condition 

[(1  -  o2z2)  -  £]  =  0  or  z  =  ±(1  -  i)1,2/a 
and  the  excursion  is  twice  this  as  expected. 


3.3  Implementation 

In  order  to  determine  the  sound  speed  with  these  methods,  we  require  a  set  of  measurements 
of  the  horizontal  wavenumber  in  the  waveguide  and  the  sound  speed  for  at  least  one  depth  ( 
this  last  permits  conversion  of  the  index  of  refraction  profile  to  a  sound  speed  profile).  For 
inversion  method  2,  we  assume  that  the  sound  speed  at  the  surface  and  both  the  depth  an 
value  of  the  minimum  sound  speed  in  the  profile  are  known.  From  the  calculations  of  the  WKB 
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phase  integral  values  for  the  various  cases,  it  is  clear  that  the  mode  numbers  of  the  measured 
eigenvalues  must  be  correctly  identified,  and  that  the  boundary  interactions  experienced  by 
each  mode  must  be  correctly  identified.  From  the  eigenvalues  and  associated  WKB  phase 
integral  values,  the  variables  {  =  (kTjl/k o)2  and  F(£)  are  computed.  Une  additional  point, 
corresponding  to  horizontally  propagating  energy,  is  added  at  £  =  1  and  F(£  =  0)  =  0. 

The  first  step  in  either  inversion  method  is  to  fit  the  data  using  either  a  least  squares  basis 
spline  or  a  cubic  smoothing  spline.  Either  of  these  two  spline  methods  allows  for  noise  in  the 
data.  The  spline  routines  essentially  provide  a  piecewise  polynomial  that  provides  a  best  fit 
to  the  data  in  a  sense  which  depends  on  the  type  of  spline  routine  used.  The  use  of  splines  is 
convenient  in  that  the  integrations  and  differentiations  required  for  the  inversion  are  easily 
performed  once  the  coefficients  for  the  spline  representaion  have  been  computed.  All  the 
computer  routines  used  for  fitting  the  data  with  the  spline,  differentiation,  and  integration 
are  in  the  IMSL,  Inc.  MATH/LIBRARY  [1]  . 

Before  discussing  the  differences  in  the  two  spline  fitting  methods,  it  is  worthwhile  to 
define  the  basis  spline  which  first  requires  a  definition  of  the  divided  difference  (for  a  com¬ 
plete  discussion  of  splines  see  reference  [14]).  For  a  function  g  which  is  given  at  the  points 
x i, . . .  the  kth  divided  difference  of  g  at  the  points  x,, . . .  ,x,+i  is  the  leading  coefficient 
of  the  polynomial  of  order  k  +  1  which  agrees  with  g  at  the  points  x,, . . . ,  xI+t;  it  is  denoted 
by  [x,, . . . ,  x1+*]#.  The  agreement  between  the  two  functions  referred  to  here  means  that  if  a 
point  occurs  m  times  in  the  sequence  x  ,  then  the  two  functions  and  their  derivatives  agree 
m-fold  in  that  the  two  functions  and  the  derivatives  up  to  the  (m-  1  )th  derivative  are  equal. 
For  a  sequence  (x0,xi,. . .  ,xyv),  if  we  let 

9jv(x)  =  (x  -  x0) . . .  (x  -  xn) 

and 

®  Ar(®»)  =  (*i  —  *0)  •  •  •  (Xj  —  X\—  i  )(X,  —  X,’+i  )  . .  .  (x,'  —  Xtf) 

(i.e.  the  (x  -  xt)  factor  is  removed  before  evaluating  the  product  at  the  point  x,),  then  the 
divided  difference  is  [4] 

[xo, . . .  ,xjv]<?  =  ^  /_  \  (3.36) 

>=o  VN\X)) 
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The  basis  spline  or  B-spline  is  defined  as  follows  [14]:  for  a  nondecreasing  sequence  t  =  (t,), 
the  ith  (normalized)  B-spline  of  order  k  for  the  knot  sequence  t  is  denoted  by  and  is 

defined  by  the  rule 

bi,k,t  =  (^i+fc  —  ^«)[^<> •  •  •  • 

The  truncation  function  ( x  -  t)+  is  defined  as  max{0,  x  -  f}  and  the  ‘notation  is  used  to 
indicate  that  i  is  fixed  and  ( t  -  x)*_1  is  to  be  considered  as  a  function  of  t  alone.  Note  that 
Bi,k,t(x)  o  for  X  g  [t, ,*,+*].  A  spline  function  of  order  k  with  knot  sequence  t  is  any  linear 
combination  of  B-splines  of  order  k  for  the  knot  sequence  t. 

The  least  squares  spline  routine  calculates  a  weighted  discrete  Z-2  approximation  to  the 
given  data  (x,, /,)  for  i  =  1,2,..., TV,  (it  finds  the  coefficients  a3)  to  minimize  the  weighted 
square  error  between  the  data  and  the  spline  i.e. 

£  I  /.  -  Haj£j(x.)  !2  m- 

>=i 

The  number  of  data  points  is  N,  B}  is  the  jth  spline  of  order  k,  m  is  the  number  of  coefficients 
(the  number  of  B-splines  making  up  the  polynomial  representation),  and  wt  are  the  data 
weights.  The  order  k  is  the  order  of  the  polynomial  pieces  (  a  polynomial  of  order  k  is  a 
polynomial  of  degree  k  -  1;  for  a  cubic  polynomial  k  =  4).  In  general,  the  weights  for  our 
problem  have  been  selected  such  that  the  measured  data  are  equally  weighted,  and  the  added 
point  at  (1,0)  is  very  heavily  weighted  compared  to  the  measurements.  These  weights  can 
be  tailored  to  fit  the  confidence  in  the  measurements. 

The  smoothing  spline  is  a  natural  cubic  spline  approximation  with  knots  at  the  data 
abscissas  where  the  term  ‘natural’  refers  to  the  end  condition  imposed.  In  addition  to  the 
constraints  imposed  by  the  required  agreement  between  the  function  and  the  spline  at  the 
knots,  some  conditions  or  constraints  must  be  given  at  the  endpoints  of  the  interval  on  which 
the  spline  is  being  calculated.  In  the  case  of  the  natural  spline  of  order  k  =  2  and  N  data 
points,  the  condition  specified  is 

i )  =  f(l){xN)  =  0;  j  —  m . k  -  2 

which  in  the  case  of  a  cubic  spline  (m  =  2)  becomes 

/(2)(xi)  =  f{2)(xs)  =  0. 
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In  our  problem,  the  data  error  variance  is  unknown,  and  the  required  amount  of  smoothing 
is  determined  in  the  IMSL  smoothing  spline  routine  by  using  the  generalized  cross  validation 
method  of  Craven  and  Wahba  [13].  Modelling  the  data  as 

y(<)  =  p(<)  +  «(0;  <€[o,i], 

the  problem  is  to  construct  g(t)  from  the  discrete  data  points  The  function  g(t)  is  assumed 
to  be  in  W^,  where 

=  {g  :  gv  abs.  cont.,  v  —  0, 1, . . .  ,m  -  1,  g^  G  /^[0, 1]}. 

The  estimate  of  g  is  gn<\,  the  solution  to  the  problem:  Find/  G  W ^  to  minimize 

j|E  </(<-> +  (/<->(«))’*«• 

The  solution  gn<\  is  the  smoothing  spline  of  order  2m- 1  [38]  with  the  parameter  A  controlling 
the  tradeoff  between  the  smoothness  of  the  solution  measured  by 

fQ  [/(m)(«)]2d« 

and  the  infidelity  to  the  data  measured  by 

;E  (/«-)-»)’• 

j- 1 

The  generalized  cross-validation  estimate  for  the  minimization  of  the  average  square  error  is 
defined  by  [13] 

V( A)  =  \  ||  (/  -  -4(A))  y  ||2  /  |~Trace(J  -  -4(A))] *  ,  (3.37) 

where  y  =  (yi,. . . ,  ynY  and  -4(A)  is  the  n  x  n  matrix  satisfying 

(9n,x(h ),...,  gn,\{tn)y  =  A(A)y. 

Once  the  initial  spline  fit  to  the  data  has  been  done  using  either  of  the  two  methods,  the 
paths  of  the  two  inversion  methods  diverge.  The  important  points  related  to  splines  for  this 
work  are:  (1)  the  spline  routines  fit  piecewise  polynomials  to  the  data  (the  polynomial  pieces 
are  generally  cubic,  although  the  least  squares  routine  allows  the  order  of  the  polynomial 
pieces  to  be  specified);  (2)  the  spline  routines  used  allow  foT  the  presence  of  noise  in  the  data; 
and  (3)  the  integrations  and  differentiations  are  based  on  the  polynomial  coefficients  from 
the  spline  fit  to  the  data. 
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3.3.1  Implementation  of  inversion  method  1 


The  function  /f(£)  is  generated  by  taking  the  derivative  of  the  spline  representation  of  F(£) 
at  the  data  points  £,  of  the  input  data.  A  spline  is  then  fitted  to  //(£).  Since  noise  in  the 
original  data  may  be  amplified  in  the  processs  of  computing  the  derivatives,  a  smoothing  or 
least  squares  spline  is  also  used  at  this  stage  of  the  computations.  The  integral  that  must 
now  be  evaluated  is 


=  _  2  ['  H(Qdt 

■K  JK  (£  -  K 


(£  - 


(3.38) 


which  has  a  singularity  at  k.  To  avoid  problems  with  the  singularity  during  integration,  the 
change  of  variables  ([4],  [43])  tz2  =  (£  -  n)  is  made  and  the  integral  to  be  evaluated  becomes 


4  /( 

=  —  /  H(k  +  u2)du. 

*  Jo 


(3.39) 


This  is  evaluated  at  points  over  the  range  of  k  (which  coincides  with  the  range  of  £)  using 
the  spline  and  integration  routines.  The  one  required  measurement  of  c(z)  is  then  used  to 
convert  the  points  of  (z,  71(2))  to  a  sound  speed  profile  c(z). 


3.3.2  Implementation  of  inversion  method  2 

The  integral  which  must  be  evaluated  to  obtain  the  inclusion  is 

2  rK  F(S') 

W*  (3'40) 

which  is  of  the  same  form  as  equation  ( 3.38),  and  the  same  change  of  variables  is  used  to 
eliminate  the  singularity.  The  integration  is  again  done  using  the  same  spline  routines.  The 
resulting  points  of  the  function  /(£)  are  fitted  with  a  spline  and  differentiated  with  respect 
to  £  to  produce  the  desired  result,  -A(£). 


3.4  Application  of  the  Inversion  Methods  to  the  Analytic 
Profiles 

The  two  analytically  solvable  profiles  introduced  earlier  provide  excellent  examples  with  which 
to  test  the  inversion  methods  since  the  results  of  the  inversion  calculations  can  be  compared 
to  known  values  during  the  intermediate  stages  of  the  computations.  Expressions  for  the 
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eigenvalues  of  these  sound  speed  profiles  are  well  known,  and  the  eigenvalues  can  be  computed 
without  the  use  of  a  propagation  modelling  program.  In  the  case  of  the  n2(z)  linear  profile,  we 
take  the  ocean  as  a  half  space  with  a  pressure  release  surface,  and  use  a  =  2.28  x  10-5  m-1 
(this  gives  a  sound  speed  gradient  of  0.017sec-1  at  a  depth  of  1000  meters  corresponding 
to  sound  speed  profile  where  pressure  is  the  only  factor  affecting  the  sound  speed  i.e.  an 
isovelocity,  isosaline  ocean).  For  the  parabolic  profile  the  ocean  is  treated  as  an  infinite 
medium,  and  the  coefficient  a  in  the  equation  for  n2(z)  is  1.257  x  10-4m“2.  In  each  case 
the  50  eigenvalues  used  have  been  calculated  for  a  frequency  of  220  Hz  and  a  1500  meter  per 
second  reference  sound  speed.  Unless  otherwise  specified,  all  the  calculations  have  been  done 
using  the  least  squares  basis  splines. 

3.4.1  Method  1  applied  to  the  n2(z)  linear  profile 

The  data  points  for  the  function  F(£)  are  created  by  dividing  the  input  values  of  the  horizontal 
wavenumber  by  ko  =  w/co,  and  squaring  the  result  to  give  £.  The  eigenvalues  for  this  example 
were  calculated  using  the  expression  [8] 

*./**  (3.41) 

where  yn  are  the  Airy  function  zeros,  H  =  (afc£)-1^3  (equation(2.28)),  and  a  is  the  coefficient 
in  the  equation  for  the  index  of  refraction.  The  asymptotic  expressions  [2]  were  used  for  the 
Airy  function  zeros;  for  the  particular  parameters  used  in  thes  examples  this  approximation 
did  not  introduce  a  significant  error.  The  corresponding  value  of  F(f )  is  obtained  by  dividing 
the  assumed  value  of  the  WKB  phase  integral  by  ko.  For  this  profile,  all  modes  interact  with 
the  surface  and  have  one  turning  point,  and  the  WKB  phase  integral  is  given  by  equation 
(3.4) 

/  (*-  01,2dz  -r{n~  1/4);  n=  1,2,3,.... 

Jo 

The  point  at  (1,0)  is  added  to  the  data  set,  and  the  spline  coefficients  calculated.  Figure 
(  3-2)  shows  the  function  F(£)  for  this  particular  case.  The  background  curve  represents 
F(f)  computed  using  the  analytic  expression,  and  the  discrete  points  are  the  values  from 
the  inversion  calculations  showing  that  there  is  good  agreement  between  the  actual  value  of 
the  phase  integral  and  the  assumed  value.  Using  the  spline  coefficients,  the  first  derivative 
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Figure  3-2:  F({)  for  the  n2(z)  linear  profile.  The  discrete  points  are  the  input  data  to  the 
inversion  program  generated  from  the  horizontal  wavenumbers  and  the  associated  values  of 
the  WKB  phase  integral.  The  curve  is  F(()  computed  using  the  analytic  expression. 

of  F(()  with  respect  to  £  is  computed  giving  H({)  as  shown  in  figure  ( 3-3).  In  calculating 
H(£),  a  noticeable  error  occurs  near  £  =  1.0  due  to  the  end  conditions  in  the  spline  fit  for 
F(£)  which  do  not  coincide  exactly  to  the  behaviour  of  the  function. 

A  spline  is  fitted  to  the  function  H(£)  represented  by  the  the  resulting  discrete  values  of 
H(Z),  and,  using  the  change  of  variables  u2  =  (£  -  k),  the  integral 

z  —  —  /  H(k  +  u2)du 

5T  Jo 

is  evaluated  by  varying  n  through  the  range  for  which  input  data  is  available.  In  essence,  for 
each  chosen  value  of  k,  the  depth  at  which  a  ray  with  the  local  propagation  angle  defined 
by  k  is  found.  As  a  matter  of  convenience  the  range  of  k  or  {  has  been  equally  divided  into 
the  number  of  equal  segments  corresponding  to  the  desired  number  of  depth  values  to  be 
computed.  Figure  ( 3-4)  shows  the  results  of  the  inversion  for  the  n2(z)  linear  profile  using 
50  modes  calculated  at  a  frequency  of  220  Hz.  The  convention  in  presenting  sound  speed 
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Figure  3-3:  H(£)  for  the  n2(z)  linear  profile.  The  discrete  points  are  the  values  of  //(£) 
computed  by  the  inversion  program.  The  curve  is  the  function  calculated  using  the  analytic 
exp-ession. 

profiles  will  be  to  show  inversion  outputs  as  discrete  points  with  the  true  profile  plotted  as 
a  solid  curve.  The  sound  speed  error  as  a  function  of  depth  for  this  case  is  shown  in  figure 
(3-5).  Note  that  the  largest  errors  occur  near  the  surface  where  the  spacing  between  the  data 
points  along  the  £  axis  is  greatest,  and  where  the  extrapolation  from  the  measured  data  to 
the  added  point  placed  at  (  =  0  is  made.  This  region  also  corresponds  to  the  region  with  the 
greatest  error  in  the  calculated  values  of  B((). 

3.4.2  Method  2  applied  to  the  n2(z)  linear  profile 

The  initial  data  preparation  is  the  same  for  both  inversion  methods  since  this  is  a  monotonic 
profile.  Once  the  spline  has  been  computed  for  the  function  F(£),  the  inclusion  is  computed 
from  equation  ( 3.40)  using  the  change  of  variables  u2  =  -  (  to  avoid  problems  with  the 
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singularity  at  ('  =  £.  The  integral  to  be  evaluated  becomes 

4 

7(0=7/  F{u2  +  t)du 

T  Jo 

After  computing  a  spline  fit  to  the  points  of  /(£)  generated  by  the  integration,  we  take  minus 
the  derivative  of  /(f)  with  respect  to  f  to  arrive  at  the  excursion  as  shown  in  figure  ( 3-6). 
The  error,  plotted  as  sound  speed  error  as  a  function  of  depth,  is  shown  in  figure  ( 3-7). 

3.4.3  Method  1  applied  to  the  parabolic  profile 

The  values  of  f  for  the  given  input  data  are  calculated  in  the  same  way  as  for  the  previous 
case;  however,  the  values  of  F( f )  are  calculated  using  equation  ( 3.6) 

k0  r  (*-01/2  dz  =  it  (n-  1/2);  n  =  1,2,3 . 

Jz,7 

since  there  are  two  turning  points  and  no  boundary  interactions  for  this  profile.  Because 
the  turning  depths  are  equidistant  from  the  channel  axis,  F(f)  is  divided  by  two,  and  the 
inversion  proceeds  as  before.  The  resulting  sound  speed  profile  is  for  the  region  below  the 
channel  axis,  and  it  is  simply  a  matter  of  using  the  profile  symmetry  to  generate  the  other 
half  of  the  profile.  Figure  ( 3-8)  shows  the  results  of  applying  inversion  method  1  to  the 
parabolic  profile,  and  figure  ( 3-9)  is  plot  of  the  inversion  error  presented  as  excursion  error 
as  a  function  of  depth.  Excursion  error  as  a  function  of  depth  is  a  clearer  way  to  evaluate 
the  inversion  results  for  two  branch  profiles  recovered  using  method  2,  and  for  consistency 
the  errors  for  inversion  method  1  are  presented  in  the  same  manner.  The  excursion  error  in 
the  case  of  a  symmetric  profile  recovered  using  inversion  method  1  is  twice  the  error  of  the 
output  for  the  single  branch  of  the  profile. 

3.4.4  Method  2  applied  to  the  parabolic  profile 

The  parabolic  profile  involves  modal  interactions  with  two  turning  points,  and  inversion 
method  2  gives  only  excursion  as  a  function  of  sound  speed.  Given  prior  knowledge  that 
the  profile  is  symmetric  this  may,  if  desired,  be  converted  to  a  sound  speed  profile  using 
symmetry  as  in  inversion  method  1.  Figure  ( 3-10)  is  the  excursion  for  the  parabolic  profile 
computed  using  inversion  method  2,  and  figure  ( 3-11)  plots  excursion  error  as  a  function  of 
sound  speed  for  this  inversion. 
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Excursion  error  (m) 


3.5  Application  of  the  inversion  methods  to  a  bilinear  profile 


The  sound  speed  profile  in  the  ocean  is  determined  by  the  temperature,  pressure,  and  salinity 
which  in  general  leads  to  sound  speed  profiles  more  complicated  than  the  ideal  profiles  we  have 
considered.  A  nominal  deep  ocean  sound  speed  profile  consists  of  a  region  near  the  surface, 
the  thermocline,  where  the  temperature  effects  dominate  and  the  sound  speed  decreases  from 
its  surface  value  to  a  minimum  at  the  axis  of  the  deep  sound  channel.  Below  the  depth  of 
minimum  sound  speed,  the  increase  of  pressure  with  depth  is  the  dominant  influence,  and 
the  sound  speed  increases  with  depth.  Although  the  sound  speed  structure  can  be  quite 
complicated,  we  will  use  a  simple  bilinear  model  with  a  channel  axis  at  686  meters  (figure 
( 3-12)  to  further  illustrate  the  inversion  methods. 

The  modal  eigenvalues  for  this  model  were  computed  for  a  frequency  of  220  Hz  using  the 
SACLANTCEN  normal-mode  acoustic  propagation  model  program  (SNAP)[27).  Initially 
the  data  set  was  truncated  before  applying  inversion  method  2  in  order  to  use  only  modes 
whose  upper  turning  point  is  deeper  than  100  meters  (to  avoid  effects  of  interactions  with  the 
surface).  The  resulting  excursion  as  a  function  of  sound  speed  is  shown  in  figure  (  3-13),  and 
the  excursion  error  is  shown  in  figure  (3-14).  The  greatest  error  occurs  near  the  channel  axis 
where  the  profile  smoothness  assumptions  are  not  met,  and  where  the  extrapolation  from  the 
measured  data  to  the  assumed  point  at  £  =  1.0  is  made. 

Next  we  use  a  data  set  containing  the  eigenvalues  for  normal  modes  with  a  lower  turning 
depth  shallower  than  3200  meters.  The  initial  portion  of  this  data  set  consists  of  eigenvalues 
for  modes  with  two  turning  points  while  the  eigenvalues  at  the  end  of  the  data  set  are  from 
modes  that  interact  with  the  surface  and  have  a  lower  turning  point  in  the  water  column. 
We  assume  that  the  value  of  the  sound  speed  at  the  surface  is  available  which  allows  us 
to  determine  the  value  of  the  horizontal  wavenumber  separating  the  two  types  of  modes  by 
setting  the  WKB  phase  integral  integrand  equal  to  zero  so  that 


c(z) 


-k 


Tn 


=  0. 


For  the  modes  with  two  turning  points  F({)  is  calculated  using  equation  ( 3.6) 


(3.43) 


Ftf)  =  *(»*  -  1/2 )/k0, 
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Figure  3-14:  Excursion  error  for  the  bilinear  profile  in  the  region  with  two  turning  points, 
and  for  the  modes  with  only  one  turning  point  equation  ( 3.4)  gives 

F(0  =  *(»-l/4)/fco. 

Applying  inversion  method  2  gives  the  excursion  versus  sound  speed  which  in  the  portion 
of  the  profile  having  two  turning  points  is  the  difference  in  depth  between  points  of  equal 
sound  speed,  and  in  the  portion  of  the  profile  with  a  single  turning  point  is  the  depth  at 
which  the  value  of  sound  speed  occurs.  Figure  (3-15)  shows  the  excursion  versus  sound  speed 
obtained.  As  before,  large  errors  occur  near  the  channel  axis,  plus  there  are  large  errors 
at  the  depth  where  the  transition  from  two  to  one  turning  point  is  made.  This  latter  error 
appears  to  result  from  errors  in  the  spline  fit  at  the  point  where  the  form  of  F(£)  changes.  If 
we  include  a  sediment  layer  with  a  1  m/s  sound  speed  gradient  starting  at  a  depth  of  4000 
meters,  then  the  inversion  using  method  2  will  give  the  distance  between  the  two  turning 
points  of  equal  sound  speed  for  sound  speeds  in  the  range  from  1478  m/s  to  1507  m/s  and 
the  actual  depth  at  which  the  sound  speed  occurs  for  sound  speeds  greater  than  1507  m/s 
(figure  (3-16)).  Consequently  we  obtain  the  sediment  sound  speed  profile  directly  and  over 
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a  portion  of  the  water  column  we  also  obtain  the  sound  speed  profile  directly.  Figure  (3-17) 
shows  the  excursion  error  over  the  entire  profile.  The  regions  of  largest  error  occur  where  the 
form  of  either  the  sound  speed  profile  or  the  function  F({)  changes.  The  outputs  presented 
here  assume  that  in  addition  to  the  modal  eigenvalues,  the  sound  speed  at  the  surface  and 
the  minimum  sound  speed  in  the  profile  are  available.  If  the  last  item  is  not  known,  it  may  be 
possible  to  use  an  iterative  procedure  to  find  an  approximation  to  the  excursion  as  a  function 
of  sound  speed  (i.e.  find  the  minimum  sound  speed)  or  to  separate  the  portions  of  the  sound 
speed  profiles  above  and  below  the  channel  axis;  however,  we  have  not  pursued  this.  The 
final  approach  for  this  profile  is  to  treat  it  as  a  symmetric  profile  by  defining  an  ‘equivalent 
symmetric  profile’  [35]  such  that 

*e.P  =  ±^(|*+  l  +  U-l)  (3.44) 

where  z+  and  z~  are  the  distances  above  and  below  the  channel  axis  to  points  having  sound 
speed  c.  Inversion  method  1  is  applied  to  the  data  for  the  region  with  two  turning  points 
(F({)  must  be  divided  by  2  as  in  the  case  of  the  parabolic  profile)  to  give  the  portion  of  the 
equivalent  symmetric  profile  below  the  channel  axis,  and  symmetry  is  used  to  generate  the 
portion  above  the  axis.  Figure  ( 3-18)  illustrates  the  results  of  using  this  method  which,  in 
a  sense,  is  averaging  the  vertical  phase  contributions  of  the  two  portions  of  the  sound  speed 
profile. 
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Figure  3-15:  Excursion  versus  depth  for  the  bilinear  profile  using  eigenvalues  for  modes 
turning  above  3200  meters.  Below  2813  meters  this  curve  is  the  sound  speed  profile. 


Figure  3-16:  Excursion  versus  depth  for  the  bilinear  profile.  Below  2813  meters  this  curve  is 
the  sound  speed  profile.  A  frequency  of  140  Hz  was  used  for  this  inversion  to  provide  bottom 
interacting  modes. 
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Figure  3-17:  Excursion  error  as  a  function  of  sound  speed  for  the  bilinear  profile.  For  sound 
speeds  greater  than  1507  m/s  this  is  depth  error. 
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Figure  3-18:  Sound  speed  versus  depth  for  the  symmetric  profile  equivalent  to  the  bilinear 
profile. 
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3.6  Performance  of  the  inversion  methods 


For  any  inverse  method,  questions  of  practical  interest  concerning  the  solution  are: 

•  Is  the  solution  unique  ? 

•  What  is  the  resolution  of  the  inversion? 

•  How  well  does  the  inversion  method  perform  with  imperfect  or  noisy  data? 

3.6.1  Uniqueness 

The  question  of  uniqueness  for  Sturm- Liouville  problems  has  been  extensively  studied  (  see  for 
example  [29],  [32]),  and  the  requirements  for  unique  recovery  of  the  coefficients  are  well  known. 
If  we  model  the  problem  as  a  horizontally  stratified  medium  consisting  of  the  water  column, 
fluid  sediment  layers,  and  a  semi-infinite  halfspace  which  can  be  considered  as  providing  an 

impedance  boundary  condition,  then  the  problem  can  be  written  on  a  finite  interval  [0, 1]  as 

^  +  (A -«(*))*>  =  0  (3.45) 

sin  at>(0)  +  co6  at/(0)  =  0  (3.46) 

sin/?v(l)  +  cos/?t/(l)  =  0  (3.47) 

The  general  problem  is  to  find  a,/3,  and  q  from  spectral  data.  The  basic  data  set  is 
the  infinite  set  of  eigenvalues  Aj  <  As  <  ...  for  the  given  Sturm- Liouville  problem.  If  the 
potential  is  symmetric  about  s  =  1/2  and  a  =  sr  -  /?  this  is  sufficient  to  provide  for  unique 
recovery  of  a  and  q;  however,  for  a  general  profile,  an  additional  infinite  sequence  of  data  will 
be  required.  Three  different  possibilities  for  this  second  set  of  data  include  [32] 

1.  a  second  set  of  eigenvalues  Aj  <  A2  <  ...  where  a  different  boundary  condition  is 
applied  in  place  of  equation  ( 3.47),  i.e.  /?  would  be  replaced  by  ft; 

2.  the  set  of  normalization  constants  6n  =||  vn  ||2  /[®n(0)]2,  n  =  1,2,...,  when  0  <  a  <  7r 
or  6n  =||  vn  ||2  /[vj*(0)]2  when  a  =  0.  is  the  eigenfunction  for  the  eigenvalue  A„,  and 

II  Il’=  tilvn(s)?d,; 
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3.  the  set  of  ratios  vn(l)/t>„(0)  when  0  <  a, (3  <  r  or  v{,(l)/^(0)  when  a  =  /?  =  0,  or 
similar  ratios  for  a  =  0,0  <  /?  <  t  or  0  =  0,0  <  a  <  t. 

In  addition,  Hochstadt  and  Lieberman  [26]  show  that  uniqueness  is  possible  given  the  set  of 
eigenvalues  Ai  <  Xj  <  . . .  and  knowledge  of  q(s)  on  the  interval  [1/2, 1]. 

For  the  problem  considered  here,  we  have  one  partial  set  of  eigenvalues  and  knowledge  of 
one  boundary  condition  (the  pressure  release  surface).  Conseqently,  we  cannot  expect  that 
the  answer  provided  by  the  inverse  methods  will  be  unique.  From  a  practical  standpoint  this 
means  that  there  may  be  additional  information  on  a  finer  verticle  scale  than  the  frequency  of 
our  experiment  could  measure,  but  which  could  be  obtained  at  a  different  frequency;  or  there 
may  be  a  variety  of  models  possible  for  the  region  below  the  depth  of  deepest  penetration 
of  the  rays  equivalent  to  the  modes  used  for  the  inversion  which  would  give  the  same  finite 
set  of  eigenvalues  as  we  have  measured.  It  should  be  noted  that,  if  horizontal  wavenumbers 
for  two  (or  more)  frequencies  are  used,  the  transformation  to  the  variable  £  will  map  all  the 
eigenvalue  data  to  the  inerval  [0,1],  and  all  the  data  will  fit  on  the  same  curve  F(£)  since 
F(f)  is  normalized  This  allows  us  to  combine  data  obtained  using  two  frequencies. 

3.6.2  Resolution 

One  of  the  advantages  of  perturbative  inversion  methods  is  that  the  linear  inverse  theory 
used  allows  a  quantitative  estimate  of  the  resolution  provided  by  the  data  (for  a  discussion  of 
linear  inverse  theory  see  [33]).  With  these  nonlinear  inversion  methods  we  can  only  provide 
a  qualitative  description  of  the  depth  resolution.  Inversion  method  1  uses  the  function  H{£) 
obtained  by  interpolating  and  differentiating  the  input  data  to  find  the  depth  at  which  a 
particular  value  of  the  index  of  refraction  will  cause  a  ray  with  angle  sin-1  to  turn.  We 
can  expect  that  the  resolution  will  be  related  to  the  depth  difference  between  turning  depths 
for  successive  modes.  Better  resolution  should  be  obtained  if  the  sound  speed  profile  and 
the  frequency  are  such  that  the  turning  depths  of  adjacent  modes  are  closely  spaced.  The 
turning  depths  are  found  using  the  condition 

33  =*r"’ 

and  high  frequencies  or  large  sound  speed  gradients  give  more  closely  spaced  turning  depths. 
In  integrating  over  k  to  obtain  the  depth  for  a  given  index  of  refraction,  we  divided  the  range 
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of  k  uniformly,  and  the  output  values  from  the  inversion  will  be  closer  together  in  portions 
of  the  sound  speed  profile  with  the  highest  gradients.  This  can  be  seen  in  figure  (  3-16). 

Inversion  method  2  approaches  the  problem  in  a  similar  manner  finding  the  difference 
in  depth  between  points  on  the  profile  having  the  same  sound  speed,  and  similar  arguments 
concerning  the  spacing  of  the  mode  turning  points  pertain  to  the  depth  resolution  of  this 
method. 


3.6.3  Performance  of  the  inversion  methods  in  the  presence  of  noise 

Because  of  the  nonlinearity  of  the  problem,  it  is  difficult  to  obtain  an  analytic  expression 
giving  an  indication  of  how  the  inversions  will  perform  when  the  input  data  is  not  exact.  For 
inversion  method  1,  we  obtained  (equation  (3.18))  the  Abel  integral  equation 

mo  _  H(n  _  _  i  r  dz 

di  2  Jo  [*(*)- flV2 

and  a  similar  Abel  equation  was  obtained  in  method  2  where,  rather  than  inverting  the 
equation,  we  differentiated  with  respect  to  £  to  obtain  the  desired  result,  equation  (  3.29) 
(the  excursion).  We  have  consideied  that  the  values  of  F(£)  are  exact,  and,  while  some  noise 
will  be  introduced  in  differentiation  due  to  the  inexactness  of  the  values  {  associated  with 
F(£),  the  noise  enters  primarily  in  the  kernel  of  the  integral  equation.  Thus  we  will  approach 
the  question  of  noise  effects  by  numerically  testing  the  inversions  using  the  n2(z)  linear  profile 
as  a  test  case. 

In  carrying  out  an  experiment  to  measure  the  eigenvalues  in  the  ocean  waveguide,  the 
environment  and  the  instruments  will  introduce  noise  into  the  data,  the  process  of  computing 
the  eigenvalues  will  have  some  inherent  resolution  limit  imposed  by  array  size  and  spatial 
sampling  factors,  and  the  computation  process  will  introduce  numeric  inaccuracies  such  as 
round-off  error.  Rather  than  attempt  to  quantify  all  these  factors  which  will  vary  from 
experiment  to  experiment,  we  will  assume  that  we  have  data  with  random  errors  of  known 
maximum  value.  This  noise  will  be  modelled  as  am  additive  roundoff  error  with  a  uniform 
probability  distribution  from  -5.0  x  10""  to  5.0  X  10“"  where  n  is  the  lairgest  decimal  position 
in  the  eigenvalue  affected  by  the  noise.  The  two  inversion  methods  were  run  with  varying 
error  magnitudes  using  both  the  least  squares  basis  spline  routines  and  the  cubic  smoothing 
spline  routines  relying  on  cross  validation  to  choose  the  smoothing  parameter.  Figures  (3-19) 
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25  0 


Log  max  data  error 

Figure  3-19:  RMS  error  for  inversion  method  1  using  least  squares  basis  splines  with  noisy 
data. 

through  figure  (3-26)  plot  the  root  mean  square  error  and  the  absolute  value  of  the  maximum 
inversion  errror  as  functions  of  the  logarithm  of  maximum  introduced  error.  The  results  show 
that  the  smoothing  spline  can  handle  a  greater  magnitude  of  noise,  and  that  both  inversion 
methods  can  handle  reasonable  amounts  of  noise  in  the  data.  If  a  mode  is  missed, 

i.e.  the  mode  spacing  and  noise  are  such  that  one  mode  is  missed  and  all  subsequent  modes 
are  incorrectly  numbered  by  one,  then  the  errors  in  F(£)  will  produce  errors  in  the  output. 
Figure  ( 3-27)  illustrates  the  error  produced  when  the  eigenvalue  for  mode  25  of  the  n2(z) 
linear  profile  is  deleted,  and  the  remaining  24  eigenvalues  misnumbered.  The  sudden  jump 
in  F(Z)  leads  to  a  sharp  jump  in  the  error  of  the  output  at  the  depth  corresponding  to  the 
missing  mode,  and  then  the  error  decreases  tending  to  a  constant  value. 
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Log  max  data  error 


Figure  3-22:  Absolute  value  of  the  maximum  error  for  inversion  method  2  using  least  squares 
basis  splines  with  noisy  data. 


Figure  3-23:  RMS  error  for  inversion  method  1  using  smoothing  splines  with  noisy  data. 


Figure  3-24:  Absolute  value  of  the  maximum  error  for  inversion  method  1  using  smoothing 
splines  with  noisy  data,. 


Log  max  data  error 


Figure  3-25:  RMS  error  for  inversion  method  2  using  smoothing  splines  with  noisy  data. 
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Log  max  data  error 


Figure  3-26:  Absolute  value  of  the  maximum  error  for  inversion  method  2  using  smoothing 
splines  with  noisy  data. 


Figure  3-27:  Sound  speed  error  for  the  n2(z)  linear  profile  with  mode  25  eigenvalue  removed, 
and  higher  modes  incorrectly  identified. 


Chapter  4 

Application  to  a  Shallow  Water 
Waveguide 

In  order  to  illustrate  the  application  of  this  inversion  technique  to  a  specific  problem  in 
underwater  acoustics,  we  will  consider  the  problem  of  determining  the  acoustic  properties 
of  the  sediments  in  a  shallow  water  waveguide.  The  objective  is  to  provide  a  portion  of  the 
geoacoustic  model  for  the  waveguide.  A  complete  geoacoustic  model  [25]  includes  information 
on  the  geology  and  topography  of  the  bottom,  density  profiles,  shear  wave  properties,  and 
compressional  wave  properties.  This  inverse  method  is  capable  of  providing  estimates  for 
compressions!  wave  speed  and  density  profiles,  although  we  can  only  illustrate  application  to 
finding  compressional  wave  speed  profiles. 

The  data  for  the  inversion  is  obtained  from  an  experiment  [20],  [30]  in  which  a  continuous 
wave  source  at  a  fixed  depth  is  towed  away  from  a  pair  of  moored  receivers  which  record 
sound  field  pressure  (amplitude  and  phase)  as  a  function  of  range  (see  figure  4.1).  In  actually 
carrying  out  the  experiment,  one  source  operating  at  two  different  frequencies  is  used  which, 
in  principle,  allows  determination  of  the  density  as  well  as  sound  speed.  The  pressure  field  is 
sampled  at  least  every  half  wavelength  in  range  over  an  aperture  of  several  kilometers,  and 
estimates  of  the  horizontal  wavenumbers  for  the  propagating  sound  field  are  obtained  from 
the  pressure  field  measurements  either  by  using  Prony’s  method  [15]  which  models  the  field  as 
a  sum  of  complex  exponentials,  or  by  numerically  Hankel  transforming  the  measurements  to 
obtain  the  depth  dependent  Green’s  function.  In  previous  work,  the  bottom  properties  have 
been  determined  using  the  data  from  this  type  of  experiment  [21]  by  assuming  a  geoacoustic 
model,  computing  the  theoretical  Green’s  function  for  the  model,  and  calculating  the  mean 


69 


13.9  m 


Radar  Ranging  System 


SOURCE 
Zo  =  6.1  m 
f  ,  *  140  Hz 
f  2  *220  Hz 


■Synthetic  Aperture 


C  *  1500  m/s 
p  =  1.0  g/cm3 


Anchor 


PRESSURE 

RELEASE 

SURFACE 


SURFACE  BUOY 


j.  RECEIVER  1 
Z,  *  12.5  m 


)////!)  1  f  t  1  /  /  /  //////  /  /  7  /  7  7  /  I  /  / 

HORIZONTALLY  STRATIFIED  BOTTOM 


Figure  4-1:  Experimental  configuration  for  measuring  the  eigenvalues  in  a  shallow  water 
waveguide.  The  parameters  shown  are  for  an  experiment  in  Nantucket  Sound,  Massachusetts 
[21] 

square  difference  between  the  two  Green’s  functions.  The  model  is  varied  and  the  process 
repeated  until  the  total  mean  square  difference  over  the  four  frequency /receiver  depth  combi¬ 
nations  is  a  minimum.  Rajan  et  al.  [37]  apply  perturbation  techniques  to  the  problem  using 
the  eigenvalues  obtained  via  the  Hankel  transform  of  the  pressure  field  as  the  input. 

In  constructing  the  geoacoustic  model  with  data  from  this  experiment,  information  on 
water  column  properties  will  be  available  from  separate  measurements  to  provide  the  sound 
speed  value  needed  to  transform  the  index  of  refraction  output  into  a  sound  speed  profile. 

4.1  Bottom  Models  and  Inversions  with  Synthetic  Data 

We  will  use  two  model  waveguides  each  having  a  one  hundred  meter  deep  isovelocity  water 
column,  a  one  hundred  meter  thick  sediment  layer,  and  a  basalt  basement  (half-space)  having 
density,  compressional  wave  speed,  and  attenuation  equal  to  nominal  values  for  basalt  (shear 


70 


properties  are  not  included  in  the  models).  Two  models  will  be  used  for  the  sediment  layer: 
a  model  with  a  terrigenous  sediment  layer  made  up  of  silty  sand,  and  model  consisting  of  fine 
sand  sediment.  The  data  on  which  both  models  are  based  is  contained  in  Hamilton’s  1980 
review  article  on  geoacoustic  models  [25].  For  both  bottom  types  the  sediment  density  has 
been  taken  as  constant  due  to  constraints  in  SNAP,  the  acoustic  modelling  program  which 
was  used  to  calculate  the  horizontal  wavenumbers  for  the  sound  field  in  the  waveguides  [27]. 

4.1.1  The  Terrigenous  Bottom  Model 

The  constituents  of  terrigenous  sediments  are  derived  from  weathering  and  erosion  of  rocks 
found  mostly  on  land  [40].  Hamilton  obtained  the  following  regression  equation  for  velocity 
as  a  function  of  the  depth  in  the  sediments  using  data  assembled  for  twenty  areas  in  upper, 
unlithified  layers  of  mostly  turbidites  in  various  oceans  [25] 

V(z )  =  1.511  +  1.301;  -  0.741;2  +  0.257;3.  (4.1) 

The  depth  ;  is  in  kilometers,  and  the  velocity  V  is  in  kilometers  per  second.  The  density  in 
the  sediment  layers  has  been  taken  as  p  =  1.0  g/cm?  so  that  the  eigenvalues  will  be  governed 
solely  by  the  sound  speed  profile  and  the  boundary  conditions  at  the  ocean  surface  and  the 
sediment/basement  interface.  Water  column  sound  speed  has  been  set  equal  to  the  sediment 
sound  speed  at  the  water-sediment  interface  for  the  same  reason.  Figure  ( 4-2)  is  the  sound 
speed  profile  for  the  waveguide  with  a  terrigenous  bottom  model. 

SNAP  produced  eigenvalues  for  53  propagating  modes  in  this  waveguide;  however,  only 
the  first  15  of  these  eigenvalues  were  used  in  the  inversion  in  order  to  avoid  effects  arising 
from  interactions  of  the  modes  with  the  basement  half-space  since  we  are  interested  in  testing 
our  ability  to  recover  the  sediment  sound  speed  profile.  The  highest  mode  used  turned  at  a 
depth  of  176  meters.  Tne  modes  were  treated  as  one  turning  point  modes,  and  both  inversion 
methods  were  used.  Figures  ( 4-3)  through  ( 4-6)  show  the  inversion  results  with  the  errors 
plotted  as  depth  error  versus  sound  speed.  As  expected,  the  largest  errors  occur  at  the 
water-sediment  interface  where  the  extrapolation  from  the  data  point  corresponding  to  the 
first  mode  to  the  point  at  £  =  1,  F(£)  =  0  is  made  and  where  the  profile  does  not  meet  the 
smoothness  assumptions.  We  have  also  ignored  the  possibility  of  reflections  at  the  interface 
which  may  introduce  errors  in  F(£)  for  the  low  order  modes  (those  with  large  equivalent  angles 
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of  incidence).  The  apparent  increase  in  error  at  the  lower  depths  reached  by  the  inversion 
is  believed  to  result  from  inaccuracies  in  the  spline  fit  near  the  end  of  the  data.  Note  that 
the  inversion  does  not  produce  any  points  on  the  isovelocity  portion  of  the  profile  since  there 
are  no  turning  points  in  an  isovelocity  segment  (the  one  inversion  point  shown  at  the  water 
surface  is  the  reference  sound  speed  i.  e.  the  assumed  minimum  in  the  profile).  Terrigenous 
sediments  are  predicted  to  have  a  sound  speed  at  the  water /sediment  interface  lower  than 
the  water  sound  speed  [25).  The  inversion  methods  cannot  distinguish  such  a  low  velocity 
zone  because  they  are  finding  turning  depths  which  result  from  increasing  sound  speeds.  In 
figure  (  4-7)  axe  the  results  for  model  1  with  the  water  column  sound  speed  increased  to 
1545  meters  per  second  to  give  a  low  velocity  zone.  The  reference  sound  speed  used  in  the 
inversion  was  taken  as  1545  meters  per  second  which,  because  of  the  normalization  of  the 
eigenvalues,  produced  a  value  of  £  >  1  for  the  first  eigenvalue  and,  to  produce  the  results  of 
figure  (4-7),  the  first  eigenvalue  was  discarded.  Thus  the  data  give  a  clear  indication  that  a 
grossly  inaccurate  reference  sound  speed  has  been  used.  If  the  first  eigenvalue  is  retained  and 
the  correct  reference  sound  speed  of  1511  meters  per  second  used,  then  the  result  obtained  is 
shown  in  figure  ( 4-8).  These  results  indicate  that  the  inverse  methods  will  at  least  indicate 
the  presence  of  a  low  velocity  zone  at  the  water/ sediment  interface,  and  it  may  be  possible 
to  gain  some  information  about  the  value  of  the  minimum  sound  speed  through  iterative 
adjustments  to  the  reference  sound  speed. 
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Figure  4-5:  Inversion  results  for  waveguide  model  1  using  inversion  method  2.  The  back 
ground  curve  is  the  actual  sound  speed  profile  in  the  waveguide. 
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Figure  4-7:  Inversion  results  for  a  terrigenous  waveguide  with  a  low  velocity  zone  at  the 
water/sediment  interface.  The  eigenvalue  of  the  lowest  mode  was  discarded  since  it  produced 
a  value  of  (  >  1. 


4.1.2  The  Fine  Sand  Bottom 

The  fine  sand  bottom  has  a  large  discontinuity  in  the  compressional  wave  speed  at  the 
water/sediment  interface  and  large  sound  speed  gradients  in  the  top  portion  of  the  sediment 
layer.  A  regression  equation  given  by  Hamilton  for  the  sound  speed  in  the  first  twenty  meters 
of  a  fine  sand  bottom  is 

V(z)  =  I8O6.O20015.  (4.2) 

The  depth  z  is  in  meters  and  the  resulting  sound  speed  V  is  in  meters  per  second.  For 
depths  greater  than  twenty  meters,  the  profile  has  been  extended  as  a  straight  line  with  slope 
1.4168  s-1  which  is  the  sound  speed  gradient  of  equation  (4.2)  at  a  depth  of  twenty  meters. 
Figure  ( 4-9)  is  the  sound  speed  profile  for  the  fine  sand  bottom  model. 

For  this  model  SNAP  generated  eigenvalues  for  47  propagating  modes  of  which  23  were 
used  for  the  inversions  (the  turning  depth  of  the  highest  mode  used  was  174  meters).  Both 
inversion  methods  were  again  used  for  the  inversions,  and,  in  addition,  both  the  least  squares 
and  the  smoothing  spline  routines  were  used.  The  results  are  shown  in  figures  (4-10)  through 
(4-17).  An  attempt  was  made  to  treat  the  modes  with  equivalent  incidence  angles  greater 
than  the  critical  angle  as  interacting  with  a  critically  reflecting  bottom  which  was  modeled 
as  a  1727  meter  per  second  half-space  (i.e.  they  were  treated  as  case  3  modes);  however, 
this  made  the  agreement  between  the  inversion  results  and  the  actual  profile  worse  than  for 
previous  inversions  which  treated  all  the  modes  as  simple  one  turning  point  modes  (figure 
(4-18)). 

Finally,  we  used  all  the  eigenvalues  in  the  inversion  process  for  this  model  (figure  (4-19)). 
The  inversion  output  sound  speed  approaches  the  basement  sound  speed  value  while  the 
depth  value  does  not  depart  significantly  from  the  sediment /basement  interface  depth. 
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BASEMENT:  p  =2  3  g/cm  3 
c  =  4000  0  m/s 
a  =  0  045  dB/X 
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Figure  4-9:  The  shallow  water  waveguide  model  with  a  fine  sand  bottom 
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Figure  4-10:  Inversion  results  for  waveguide 
splines). 
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Figure  4-11:  Depth  error  for  inversion  method  1  (least  squares  splines)  applied  to  waveguide 
model  2. 


Figure  4-12:  Inversion  results  for  waveguide  model  2  using  inversion  method  1  (smoothing 
splines). 
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Figure  4-14:  Inversion  results  for  waveguide  model  2  using  inversion  method  2  (least  squares 
splines). 


Figure  4-15:  Depth  error  for  inversion  methc 
model  2. 
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Figure  4-18:  Inversion  results  for  waveguide  model  2  using  inversion  method  1  (least  squares 
splines)  treating  the  low  order  modes  as  reflecting  from  a  half-space  with  a  1727  meter  per 
second  sound  speed. 
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Figure  4-19:  Inversion  results  for  waveguide  model  2  using  inversion  method  1  (least  squares 
splines)  using  the  eigenvalues  of  all  propagating  modes. 
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4.2  Experimental  Results 


An  experiment  was  performed  in  the  Gulf  of  Mexico  near  Corpus  Christi,  Texas  during 
September,  1985  using  a  configuration  similar  to  that  of  Figure  (4.1).  For  the  Corpus  Christi 
site  the  water  depth  was  30.3  meters,  source  depth  was  22.86  meters,  and  the  two  receivers 
were  1.5  meters  and  15.0  meters  above  the  bottom  (G.  V.  Frisk  pers.  com.  1988)  Geologic 
data  [31]  for  the  region  predicts  a  17.5  meter  thick  layer  of  silty  clay  with  an  estimated  sound 
speed  at  the  water /sediment  interface  of  1517.6  meters  per  second  increasing  to  1531.5  meters 
per  second  at  the  layer  bottom.  Since  the  sound  speed  in  the  water  column  was  1545  meters 
per  second,  this  gives  a  low  velocity  zone.  Beneath  the  silty  clay  layer,  a  layer  of  very  fine 
sand  was  predicted  to  extend  to  45  meters  with  a  sound  speed  of  1737.1  meters  per  second  at 
the  interface  between  the  layers  and  increasing  to  1763.0  meters  per  second  at  the  bottom  of 
the  layer.  In  the  absence  of  better  information,  parameters  for  this  segment  were  extended  to 
a  depth  of  100  meters.  It  should  be  emphasized  that  these  values  are  based  on  a  study  of  the 
geology  of  the  region,  and  the  sound  speed  predictions  are  based  on  the  type  of  sediments  and 
the  estimates  of  Hamilton’s  work  for  such  sediment  types  [25].  Figure  (4-20)  is  the  predicted 
sound  speed  profile  in  the  sediments.  The  eigenvalues  were  obtained  by  S.  D.  Rajan  (pers. 
com.  1986)  from  the  depth  dependent  Green’s  function  of  the  measured  sound  field  via  the 
numerical  Hankel  transform,  and  are  presented  in  Table  (4.1).  Water  column  sound  speed  at 

Frequency  (Hz)  Mode  Number  Eigenvalue  (m~1) 


50.0 

1 

0.1972886 

2 

0.182872 

140.05 

1 

0.5678598 

2 

0.5572134 

3 

0.5440619 

4 

0.5252742 

5 

0.5102440 

Table  4.1:  Experimentally  measured  eigenvalues  for  the  Corpus  Christi  site. 

the  experiment  site  was  1545  meters  per  second.  Since  a  total  of  only  seven  eigenvalues  at  the 
two  frequencies  is  available,  we  combined  the  data  to  form  a  single  data  set  which  is  possible 
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Figure  4-20:  Sound  speed  profile  for  the  Corpus  Christi  area  based  on  geologic  information 
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Figure  4-21:  Corpus  Christi  area  sediment  sound  speed  profile  obtained  using  all  seven  eigen¬ 
values  at  50  Hz  and  140.05  Hz. 

since  F(£)  as  we  have  defined  it  is  frequency  independent.  The  results  of  the  inversion  using 
method  1  with  an  inversion  program  using  the  least  squares  spline  routines  of  the  PORT 
Library  [18]  are  shown  in  figure  (4-21).  (The  PORT  Library  program  was  used  because  it 
seemed  to  handle  the  small  data  set  better  than  the  program  using  the  IMSL  routines.)  The 
large  oscillation  at  the  bottom  of  the  inversion  results  suggets  a  data  error  at  the  end  of  the 
data  set  (smallest  £).  After  £  is  calculated,  the  second  50  Hz  eigenvalue  is  the  last  data  point, 
and  if  this  eigenvalue  is  removed,  the  oscillation  is  also  removed  as  shown  in  figure  ( 4-22)) 
which  also  shows  the  results  of  a  perturbative  inversion  done  by  S.  D.  Rajan  (pers.  com 
1986),  and  the  predicted  sound  speed  profile. 
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Figure  4-22:  Corpus  Christi  sediment  sound  speed  profile  obtained  with  inversion  method 
1  after  removing  the  second  50  Hz  eigenvalue,  the  profile  obtained  perturbatively  by  Rajan 
(solid  curve),  and  predicted  profile  (dotted  curve). 


Chapter  5 

Conclusions 

5.1  Summary 

In  this  thesis  we  have  shown  that  using  a  WKB  inversion  technique  based  on  the  WKB 
phase  integral  equations  and  using  modal  eigenvalues  as  input  data  provides  a  method  for 
obtaining  information  about  index  of  refraction  profiles  in  the  ocean  waveguide.Information 
can  be  obtained  for  the  water  column  and  for  the  sediment  layers.  Two  approaches  for 
transforming  the  WKB  phase  integral  to  a  relation  between  depth  and  the  index  of  refraction 
were  developed.  Both  methods  assume  that  the  value  of  the  WKB  phase  integral  can  be 
accurately  estimated  based  on  the  number  of  turning  points  for  a  mode  and/or  the  types  of 
boundary  interactions  the  mode  experiences,  and  both  methods  rely  on  knowledge  of  at  least 
the  speed  of  sound  at  the  surface. 

Inversion  method  1  differentiates  the  WKB  phase  integral  equation  to  obtain  an  Abel 
integral  equation  which  is  then  inverted  to  provide  an  equation  for  depth  [43]  and  is  only 
usable  for  monotonic  or  symmetric  sound  speed  profiles.  In  the  case  of  a  symmetric  profile 
the  knowledge  that  the  phase  contributions  above  and  below  the  axis  are  equal  allow  the 
original  WKB  phase  integral  to  be  divided  into  two  equal  portions  which  essentially  puts  the 
problem  in  monotonic  form.  With  a  value  of  sound  speed  at  one  depth,  the  monotonicity 
or  symmetry  allows  conversion  of  the  depth/index  of  refraction  relation  to  a  sound  speed 
profile.  For  a  general  asymmetric  profile,  it  is  possible  to  to  treat  the  data  as  though  it  was 
obtained  from  a  symmetric  profile  and  to  obtain  an  ‘equivalent  symmetric  profile’  [35].  The 
utility  of  this  approach  remains  to  be  seen. 

Inversion  method  2  uses  fractional  calculus  methods  [39]  to  find  the  inclusion  (equation 
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(  3.28))  from  the  WKB  phase  integral.  The  inclusion  is  differentiated  to  obtain  the  depth 
difference  between  two  points  with  the  same  index  of  refraction  (the  excursion).  For  a  mono¬ 
tonic  profile  this  gives  the  depth  directly  and  for  a  symmetric  profile  the  phase  integral  can 
be  split  into  two  equal  parts  as  it  was  for  inversion  method  1  allowing  determination  of  the 
full  profile.  For  an  asymmetric  profile  this  method  gives  different  answers  in  the  portions 
of  the  profile  with  one  and  two  turning  points.  In  the  portion  with  two  turning  points,  the 
inversion  result  is  the  distance  between  two  points  of  equal  index  of  refraction,  and  in  the 
single  turning  point  case  the  result  is  the  index  of  refraction.  This  assumes  that  we  have  at 
our  disposal  the  sound  speed  at  the  surface  (to  allow  the  modes  to  be  separated  into  the  one 
and  two  turning  point  groups)  and  the  minimum  sound  speed  in  the  profile.  The  physical 
interpretation  of  the  variable  £  as  sin29  which  constrains  it  to  the  range  [0, 1],  suggests  that 
it  may  also  be  possible  to  approximate  the  minimum  sound  speed  value  based  on  the  range 
of  the  normalized  eigenvalues.  Since  a  known  portion  of  the  inversion  result  for  a  profile 
with  both  one  and  two  turning  point  modes  is  the  index  of  refraction  as  a  function  of  depth, 
it  may  be  possible  to  use  an  iterative  procedure  to  arrive  at  a  reasonable  excursion  versus 
sound  speed  curve  and,  perhaps,  to  separate  the  two  branches  of  the  sound  velocity  profile 
in  the  two  turning  point  section,  that  is  to  find  a  good  approximation  to  the  entire  sound 
speed  profile. 

Both  methods  have  been  implemented  using  splines  for  differentiation  and  integration. 
The  techniques  were  applied  to  three  idealized  profiles  to  demonstrate  the  ability  of  the 
inversion  methods  to  recover  the  profiles.  A  significant  disadvantage  of  these  methods  is  the 
present  lack  of  quantitative  statements  concerning  depth  resolution  and  performance  with 
inaccurate  data.  Performance  of  the  methods  with  inaccurate  data  was  demonstrated  with 
numerical  experiments  using  the  n2(z)  profile  after  adding  random  noise  to  the  eigenvalues. 
The  inversions  are  capable  of  giving  results  even  with  some  noise  in  the  data,  and  the  use  of 
inversion  method  1  on  a  set  of  experimental  data  was  demonstrated  in  Chapter  4.  The  results 
of  the  inversion  of  the  experimental  data  show  that  the  inversion  results  may  give  a  fairly 
clear  indication  ( large  depth  oscillations  in  the  output)  when  an  eigenvalue  is  incorrect.  The 
depth  resolution  is  related  to  the  depth  separation  between  succesive  turning  points;  however, 
we  do  not  yet  have  a  quantitative  expression  for  the  resolution. 
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In  Chapter  4,  both  methods  were  applied  to  two  models  for  a  shallow  water  waveguide  to 
illustrate  the  capability  to  provide  estimates  of  sediment  compressional  wave  speed  profiles  in 
shallow  water.  Finally,  inversion  method  1  was  applied  to  a  small  data  set  from  an  experiment 
performed  in  the  Gulf  of  Mexico.  Only  seven  eigenvalues  at  two  frequencies  were  available, 
and  it  appears  that  one  value  is  not  accurate  based  depth  oscillations  at  the  end  of  the  sound 
speed  range.  Although  our  result  is  much  different  from  the  profile  predicted  on  the  basis  of 
the  geology  in  the  area,  it  is  in  general  agreement  with  a  result  obtained  perturbatively  by 
S.  D.  Rajan  (pers.  com.  1986). 

5.2  Future  work 


A  number  of  issues  related  to  this  work  need  to  be  addressed  more  fully  in  the  future.  We  have 
largely  ignored  the  discontinuity  in  sound  speed  at  the  water/sediment  interface.  In  Chapter 
4  we  tried  treating  the  low  order  (large  angle  of  incidence)  modes  as  being  critically  reflected 
at  the  interface,  but  this  made  the  agreement  between  the  inversion  and  the  true  profile  worse 
than  ignoring  the  presence  of  the  discontinuity.  The  manner  in  which  discontinuities  in  the 
index  of  refraction  profile  should  be  treated  in  the  context  of  WKB  theory,  and  the  effect  of 
the  discontinuities  on  the  cumulative  phase  needs  further  investigation. 

In  Chapter  4  we  alluded  to  use  of  the  inversions  to  estimate  the  density  profile  as  well  as 
the  compressional  wave  speed  profile.  If  the  effects  of  density  (assumed  to  be  a  function  of 
depth  only)  are  included,  the  homogeneous  Helmholtz  equation  becomes  [8] 

V2p(r,z)  -  — ^Vp(z)  •  Vp(r,z)  +  *2(z)p(r,z)  =  0,  (5.1) 

and  separation  of  variables  gives  for  the  depth  equation 

s?™ + ■ Gpj)  zpw + (*!«  -  *0 = 0  <5-2> 

Substituting  t>(z)  =  P{z)/y/p(z)  gives  the  Schrodinger  type  equation  [36] 

•^v(z)  +  (*2(z)  +  p(z)  -  *2)  r(z)  =  0  (5.3) 


where 


/*(*)  = 


p1/2(z)  d 


d  f  1  d  ,  ; 
dz 
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The  WKB  phase  integral  is  now 


An  open  question  is  whether  we  can  obtain  a  depth  profile  for  a  new  variable  x(*)  =  n2(z)  + 
fi(z)/k0  at  two  frequencies  and  then  eliminate  n2(z)  from  the  pair  of  equations  since  the  index 
of  refraction  is  frequency  independent.  This  would  leave  an  equation  for  p(z)  which  through 
equation  (  5.4)  provides  a  differential  equation  for  p(z).  The  numerical  modelling  program 
used  in  this  work  required  constant  densities  in  the  sediments  and  the  basement  and  did  not 
provide  a  means  to  pursue  the  density  question. 

The  inverse  method  was  applied  to  one  data  set  which  had  only  six  good  eigenvalues. 
Identification  of  existing  experimental  data  from  other  experiments  obtained  with  arrays 
capable  of  resolving  a  larger  number  of  modes  (at  least  10  to  15)  would  provide  an  opportunity 
to  more  fully  test  the  inversions  with  experimental  data.  If  data  is  available  at  two  frequencies, 
use  of  Prony’s  method  [15]  to  find  the  eigenvalues  will  provide  an  estimate  of  the  modal 
attenution  in  adition  to  the  eigenvalues.  Using  these  inversion  methods  in  conjunction  with 
Prony’s  method  will  allow  a  geoacoustic  model  incorporating  sediment  sound  speed,  density, 
and  attenuation  to  be  developed. 
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